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PREFACE 

This report describes part of a comprehensive and continuing program of re- 
search in multispectral remote sensing of the environment from aircraft and 
satellites. The research is being carried out for NASA’s Lyndon B. Johnson Space 
Center, Houston, Texas, by the Environmental Research Institute of Michigan 
(formerly the Willow Run Laboratories, a unit of The University of Michigan's 
Institute cf Science and Technology). The basic objective of this program is to 
develop remote sensing as a practical tool for obtaining extensive environmental 
information quickly and economically. 

In recent times, many new applications of multispectral sensing have come 
into being. These include agricultural census -taking, detection of diseased plants, 
urban land studies, measurement of water depth, studies of air and water pollu- 
tion, and general assessment of land-use patterns. Yet the techniques employed 
remain limited by the resolution capability of a multispectral scanner. Techniques 
described in this report may help to overcome this limitation by enabling either 
examination of the contents of a given scanner resolution cell or, through averaging, 
faster estimation of the contents of a larger area. 

To date, our work on estimation of proportions has included: (1) extension 
of the signature concept to a mixture of objects; (2) development of a statistical 
and geometric model for sets and mixtures of signatures; (3) evaluation of compu- 
tational methods used to estimate proportions of a mixture by maximum likelihood; 
(4) creation of a computational technique for assessing the expected accuracy of 
estimation as a function of the signature set; (5) development of techniques to 
identify alien objects; (6) testing and evaluating the proportion estimation algo- 
rithms on artificial as well as actual multispectral scanner data; (7) examining 
the problem of establishing signatures when pure samples of the objects of in- 
terest are not available; and (8) evaluation of alternative estimators. 

The research covered in this report was performed under Contract NAS9-9784, 
Task III, and covers the period from 1 February 1973 through 31 January 1974. Dr. 
Andrew Potter has been Technical Monitor for NASA. The program was directed 
by R. R. Legault, Vice-President of the Environmental Research Institute of Michi- 
gan (ERIM); J. D. Erickson, Principal Investigator and Head of the ERIM Informa- 
tion Systems and Analysis Department; and R. F. Nalepka, Head of the ERIM Multi- 
spectral Analysis Section. The ERIM number for this report is 190100-25-T. 

The authors acknowledge the direction provided by Mr. R. R. Legault, Dr. J. D. 

Erickson, and Mr. R. F. Nalepka; the technical counsel furnished by Mr. R. J. Kauth, 
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Mr. W. A. Malila, and Dr. R. B. Crane; and the secretarial services of Miss D. 
Dickerson, Mrs. L. A. Parker, and Mrs. D. Humphrey. 


MICHIGAN 
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1 

SUMMARY 

The potential applications of remote sensing appear numerous. However, some of these 
applications are hampered by the limited spatial resolution of the sensing device. To surmount 
this difficulty, procedures have been developed for estimating proportions of objects within a 
single pixel of a multi spectral scanner. 

This report covers a third phase in the development of proportion estimation techniques. 

The first two phases included formulation of a solution to the problem, development of suitable 
algorithms for effective computation, demonstration of feasibility on more or less artificial 
data, and the pinpointing of problem areas. 

This third phase in the development of proportion estimation techniques has focused on prob- 
lems tending to limit operational usefulness. The speed of the proportion estimation program 
was increased by theoretical advances in the basic algorithm and by improvements in alien 
object detection procedures. Other attempts to increase speed included studies of averaging 
procedures for estimating proportions of objects in a group of pixels and in the use of a sim- 
plified proportion estimator. The problem of setting the alien object threshold parameter was 
studied, but an effective solution not obtained. 

In addition, proportion estimation techniques were tested on a limited amount of ERTS-1 
data where sufficient ground truth data were available. Results, though disappointing, indicate 
the need for technique modifications and further experimentation with real data sets. 
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2 

INTRODUCTION 

In recent years the stalf at ERIM has participated in the development of various tech- 
niques for using multispectral remote sensing in many applications, including agricultural 
land use measurement, rock identification, and water depth measurement. 

In conventional multispectral recognition, the total area of each ground material is mea- 
sured by identifying the material in each ground area (pixel) covered by one resolution element 
of a multispectral scanner. The total area covered by a ground material is found by adding up 
the pixels identified with that material. If almost every pixel in the ground scene contains 
just one of the possible materials, this technique provides adequate estimates of acreages. How- 
ever, if the pixel contains more than one material in substantial amounts, the pixel cannot be 
properly classified. For ERTS-1 satellite data, for example, in which each pixel covers about 
1.1 acres, the number of pixels containing more than one material may approach 30% of the total 
area for agricultural crops in the corn belt. 

The purpose of the present effort is to obtain improved area estimates of ground materials 
in these cases. We attempt to overcome the problem of boundary pixels by estimating the pro- 
portions of materials within each single pixel. 

Since its inception, this effort has consisted of a mix of theoretical model studies and tests 
with both simulated data and modest amounts of ground -truthed real data. Now that real data 
sets with adequate associated ground truth are becoming available, we will be able to utilize 
these in testing and development of mixtures procedures. The past history of the effort is 
summarized below to provide a context for this report. 

Our work on estimation of proportions was accomplished in several phases. In the first 
phase [1, 2], a mathematical model was constructed which related the multispectral signatures 
of a mixture to the signatures of component materials. This model permitted the maximum 
likelihood estimate of the proportion vector to be formulated in terms of the observed data 
point. The computational aspects of the problem required this simplification: that all of the 
covariance matrices of the signatures of the component materials be taken as equal to their 
average. Theoretical and empirical results supported the validity of this assumption. With 
this simplification, proportion estimation becomes a quadratic programming problem. Several 
existing computational methods of quadratic programming were adapted and tested on simulated 
scanner data. Results indicated that this method for proportion estimation was feasible. 

The second phase of the program [3] included investigating the problem of detecting 
alien objects — i.e,, objects in the scene not represented in the signature set. A procedure was 
devised for rejecting those pixels probably containing significant amounts of alien materials. 

In addition, aircraft scanner data were smoothed over ERTS-sized resolution elements to 
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simulate spaceborne scanner data. When proportion estimation techniques were tested on this 
data, estimates of crop acreage based on the estimated proportions were found to be better 
than estimates obtained with conventional recognition techniques. 

The third phase, covering the period of this report, was devoted to investigating problems 
in proportion estimation that limit its operational usefulness — namely computation time and 
inaccuracy stemming from the intrusion of alien objects. In addition, some preliminary trials 
were made with ERTS-1 data. Increased speed was achieved by improvement in the basic algo- 
rithm, by conversion to a more advanced computer, and by the further development of data 
averaging prior to estimating proportions. A reduction in computation time by a factor of about 
seven was attributed to improvements in the basic algorithm; the advanced computer employed 
contributed a reduction by a factor of about two. In absolute terms, it required about 20 msec 
to estimate the proportions of five materials from 12 -channel data. Averaging improved the 
speed of estimation by a factor approximately equal to the number of data points included in 
the average. {Results of theoretical analyses and simulated tests indicate that averaging may, 
under certain conditions, even improve the accuracy of proportion estimates.) With the im- 
proved proportion estimation algorithm, an ERTS frame containing about 10 pixels can be 
processed in about two hours by averaging signals in groups of 26. It would take over four hours 
to classify the pixels of the ERTS frame into seven signature categories by recognition process- 
ing using a linear decision rule. 

Also, in order to achieve increased speed, a theoretically less satisfactory but computa- 
tionally more tractable estimator of proportions was investigated. Tests on simulated data 
indicated that although computation time decreased by about 1/3, the estimate's mean-square 
error increased by about 1/3 on a pixel -by -pixel basis. A surprising result of tests with sim- 
ulated space data indicated that the new estimate could improve accuracy when estimates are 
averaged over a sufficient number of pixels (in the case of these tests, more than 27). 

In addition, during this phase, a new procedure was devised for detecting alien objects. 

This procedure depends upon setting the value of a single parameter called the alien object 
threshold. Tests with simulated data to determine a satisfactory method of setting this thresh- 
old were inconclusive. 

During the course of the year the mixtures algorithm was applied to three ERTS-1 data 
sets. Two of these — involving, respectively, lakes and rice fields — are reported elsewhere 
[4, 5] . The third case, a preliminary investigation of the use of the mixtures algorithm in a 
general agricultural situation utilizing data gathered under the CITARS (Crop Identification 
Technology Assessment for Remote Sensing) program, is reported here in Section 7. 

The results reported for lakes and for rice fields are most encouraging. These were 
achieved, however, by employing some special procedures appropriate to the situations. In the 
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case of lakes, a special detection threshold was set to keep slight indications of water in many 
pixels from accumulating to a large value. In the case of rice, mixtures were calculated only 
on those elements located on the boundaries between the rice fields. 

In the case of CITARS, the data represent a low contrast scene relative to the other two 
cases, so no special procedures were used. The results to date on this set of CITARS data are 
not impressive. It is evident that in order to develop a mixtures algorithm which will be auto- 
matic and effective when the spectral information is limited, or during seasons when contrast 
is low, additional procedures will have to be devised to complement the present algorithm. The 
only practical way to achieve this will be by utilizing large data sets with excellent ground 
truth. The past lack of such data has imposed a major constraint on the development of the 
proportion estimation technique. We are hopeful that the CITARS data will fill this need. 

The next section (Section 3) summarizes our approach to the proportion estimation prob- 
lem, describing improvements in the basic algorithm and some associated computations. Theo- 
retical and numerical results relating to proportion estimation using averaging procedures are 
presented in Section 4. In Section 5, a simplified estimator of proportions is introduced and 
compared with the standard estimator. Section 6 reviews our studies of the effect of the alien 
object threshold setting. Section 7 presents preliminary results of processing a limited amount 
of ERTS-1 data in an agricultural situation. The appendices contain proofs of theorems and 
descriptions of computer programs associated with estimation of proportions. 
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3 

APPROACH TO PROPORTION ESTIMATION 

This section summarizes ERIM’s basic approach to the proportion estimation problem. 

(See also Refs. [1, 2, 3] for more detail.) The remainder of this section is devoted to the now 
rewritten computer program (MIXMAP). Improvements made during the reporting period in 
speed of calculation are introduced in Section 3.2, and the necessary programming is detailed 
in Appendix A. Modifications of the alien object test to improve both speed and accuracy are 
described in Section 3.3. (Results of numerical trials on simulated data are reserved for Sec- 
tion 6.) Finally in Section 3.4, under the heading of signature analysis, some preliminaries to 
the basic proportion estimation calculation are described. Signature analysis provides a mea- 
sure of the expected quality of the estimates to be obtained from a given signature set. 

3.1 MODEL FOR SIGNATURES OF MIXTURES 

When the IFOV of a multispectral scanner is large with respect to the structure of the scene 
being scanned, a single resolution cell or pixel may contain more than a single object or ma- 
terial. A mathematical model has been constructed which relates the signature of a mixture of 
materials to the signatures of the component materials. Suppose the scanner has n spectral 
channels and that the signature of object class i, where 1 % i S m, is represented by the n- 
dimensional Gaussian distribution with mean A. and covariance matrix M.. Let the proportion 
of object class i be A i and let A be the vector (A^, A 2 , . . . , A m )*, where the superscript t denotes 
transpose. The signature of the mixture with proportion vector A is taken to be a Gaussian 
distribution, with mean A^ and covariance matrix given by 

\ ■ I/S ■ Ax 

M x = J\ l M. 

where A is the matrix with i column A^. These formulas constitute our model for signatures 
of mixtures of materials in terms of signatures of the individual materials. 

3.2 ESTIMATION OF PROPORTIONS 

The model for a mixture signature can be used to estimate the proportion vector A corre- 
sponding to a signal data vector from a multispectral scanner. Let y denote the n-dimensional 
data vector from the scanner. A maximum likelihood estimate of the proportion vector A [2] is 
a value of A which minimizes 

F(A) = £n|M A l + <y - A x , M^y - A^)> 
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subject to the constraints that 

FV = 1 and X 1 ^ 0 for lgigm 

Here ImI denotes the determinant of M, M _1 is its inverse, and <u, v> denotes the inner or dot 
product of the vectors u and v.* 

In general, minimizing F(X) subject to the given constraints is quite difficult. Investigations 
[2] showed that a good approximation to the minimal X could be obtained if a simplifying assump- 
tion is made. The assumption is that the average of the covariance matrices of the pure signa- 
tures can be substituted for each M . By using the simplifying assumption and applying a linear 
transformation which reduces the common covariance matrix to the identity, the problem of 
estimating X becomes one of minimizing a function G(X) of the form 

GW = lly - Ajl 2 

subject to the constraints on X. Now y represents the transformed data point, and A x the mean 
of the signature associated with the proportion vector X after the pure signature means have 
also been transformed. 

The problem of minimizing G(X) subject to the constraints on A can be viewed geometrically. 
The set of points = AX, where X is a proportion vector, is the convex hull of the A. and is 
called the signature simplex. The problem is to find a proportion vector X such that a£ is the 
point in the signature simplex closest to the data point y. 

The optimal X will be unique if the signature simplex is non -degenerate, i.e., has positive 
m - 1 dimensional volume. This is equivalent to the (n+1) -dimensional vectors ^ A*, l^j being 
linearly independent. We shall always assume this to be the case. Non-degeneracy of the 
signature simplex implies that the number of materials m in the pure signature set does not 
exceed the number of spectral channels n by more than one. 

A . 

Another estimator for X is considered in Section 5; therefore, to avoid confusion, X is 
called the ''standard" estimator. 


*Salvato [6] has pointed out the relationship of the ERIM proportion estimation model to the 
TRW mixtures model [7]. The two models are mathematically equivalent under certain condi- 
tions. 
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The problem of minimizing G(A) can be identified as a quadratic programming problem. 

In our previous work, a number of extant algorithms were adapted, the most successful based 
on the method of Theil and van de Panne. Details of the method can be found in Refs. [2, 6| . 

An improved version of this algorithm was programmed for the IBM 7094 computer during this 
contract period. The theoretical basis and programming details are given in Appendix A. This 
version is about 14 times faster than the old version which was programmed for theCDC 1604 
computer — a factor of 7 is attributable to improved methodology, and a factor of 2 is attribut- 
able to the more advanced computer. It now takes about 20 msec to estimate a mixture of 
5 materials with 12 channels of data. 


The Theil and van de Panne method adapted to the minimization of G(A) subject to the con- 
straints on X is described in Ref. [ 2, (Section 1.3)]. It depends upon projections of the data point 
y onto the linear hulls of certain subsets T of vertices A. of the transformed signature simplex. 
The linear hull of a finite set of vectors v^ . - . , v r is the set of all vectors of the form 


i=l 


.V. 


with . =1 

j=l 


When the subset of vertices contains all of the A., 1 g i g m, then the projection of y is AX where 
^ is the solution to 


*11 ■ 

• r im j 


1 

• V 

H * 

1 


g l 

r ml' 

■U 1 


A m 

= 

g 2 

1 1 . 

. 1 0 


A 


g m 


and 


( 1 ) 


r.. = < A., A.> , g. = < A., y> 
lj i j ’ & i r 

When the subset T does not contain all the vertices, the projection of y onto the linear hull of 
the vertices in T is obtained as follows. For each i for which A. is not in T, delete the i-throw 
and column of the matrix in the system of equations (1). Also delete A 1 and g. from the vectors 

and consider the remaining system. For example, if T contains all the A. but A 1 

and A.g, the linear system corresponding to T becomes: 
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The projection y onto the linear hull of the vertices in T is then given by AX where X 1 = 0 for 

A not in T, and X 1 is the value in the solution of the linear system corresponding to T for A , 
1 1 

in T. 


The new program achieves greatly increased speed by precomputing and storing the inverse 
of the coefficient matrix in the linear system corresponding to each possible subset T. Addi- 
tional efficiency is attained by utilizing recursive relationships between successive projections 
required by the Theil and van de Panne method. 

3.3 DETECTION OF ALIEN OBJECTS 

Estimating proportions of unresolved objects from a signal y is based on the assumption 
that the signal comes from a pixel which contains a mixture of materials. These materials 
are represented by known signatures that constitute the pure signature set. If the pixel should 
contain a material not represented in the signature set, significant additional error in the esti- 
mate of proportions may result. The amount of this error depends upon the proportion of these 
alien materials and the geometric relationship of their signatures to those in the pure signature 
set. Those materials occurring in a scene but not represented in the pure signature set are 
referred to as alien materials or alien objects. Procedures have been designed to reduce the 
error resulting from the presence of alien objects. These procedures take the form of thresh- 
olding tests — hence the designation "alien object threshold." 

One might attempt to avoid the alien object problem by obtaining signatures for ail materials 
present in the scene. This approachis usually impractical because of the large number of mate- 
rials present and the impossibility of obtaining definitive signatures for many of them. An alterna- 
tive is to use essentially a chi-square test as in conventional recognition processing. Some modi- 
fications are necessary when averaging procedures are also employed. 

The new mixtures program contains improved procedures for dealing with alien objects. 
These procedures can be described most easily in terms of the pure signature set and signals 
after a linear transformation has been employed. After this transformation, we assume that 
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the i-th material in the pure signature set has mean A., and its covariance matrix is the identity. 
Now given a signal (data point) y from a pixel with unknown proportions of various materials, 

A 

the estimate a of the proportion is obtained as follows. Let Z denote the point in the signature 
simplex closest to y. Then Z may be represented in the form 

A 

Z = AX 

A 

where X is a proportion vector and is taken as the estimate of proportions in the pixel repre- 
sented by the signal y. In order to apply an alien object test, we ask, "What is the probability 
that we would have observed the signal with value exceeding y if the true proportion of the pixel 

A 

was X?” Assuming Gaussian signature distributions, this amounts to a chi-square test with n 

degrees of freedom, where n is the number of spectral channels used. The level of significance 

o 

is determined by a value X Q , which is the alien object threshold. If 



then the estimate fails the chi-square test; we then say that the pixel contains significant 

amounts of alien materials and make no estimate of proportions for the pLxel in question. If 

the estimate passes the test, we accept it as the estimate of proportions of materials in the pixel 

in question. No significant theoretical results have been obtained as to where the threshold level 
2 

X Q should be set in practice, but results of empirical tests are presented in Section 6. This 
situation is analogous to the chi-square threshold setting in conventional recognition, where it was 
found that a very high threshold setting was most useful. 

In effect, the alien object test is as described for proportion estimation on a pixel- by -pixel 
basis. However, to reduce computation time, a screening test is also employed in practice. 

This screening test avoids actual computation of the proportion estimate for some fraction, of 
the cases where the estimate would fail the alien object test. The screening test is computa- 
tionally rapid and works as follows. Let Z denote the point in the signature simplex closest to 
the signal y. Set 

d 2 = lly - Zlf 2 

Let L a denote the linear hull of the signature simplex, i.e., L A is the set of points of the form 
A V where the sum of the components of rj is one. (Note that the components are not necessarily 
positive.) Let Py denote the orthogonal projection of the point y onto L A ; i.e., Py is the point in 
L a closest to y. Set 

d 2 = lly - Py II 2 
d 2 - llPy - Z II 2 
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Then 


Now let Py = At]. If the components of rj are non-negative, then set 6 = 0. If some compo- 
nent of v is negative, then Py falls outside the signature simplex and there is an extended 
simplex such that each face is parallel to its corresponding face in the signature simplex, the 
hyperplane through each face of the extended simplex is at a distance from the hyperplane 
of the corresponding signature simplex, and Py is in a face of this extended simplex. Then 


Set 


6 2 Sd d 


.2 j 2 .2 

6 =d 1 + 5 2 


Then 


« 2 S d 2 


If the alien object threshold is x^> then the screening test 

*2 ^ 2 
5 >X Q 


will eliminate points that would fail the alien object test, but will accept points that may or may 

? 2 
not fail it after d is computed. The reason for using the screening test is that 6 is easier to 

2 2 2 
compute than d , and 6 is often a good approximation to d . 


In certain averaging procedures, the screening test is used exclusively to maintain com- 
putational efficiency. It is used as follows. Let y^, . . . , y N be signals from N pixels. Of 

these N signals, let y y' be the signals which pass the alien object screening test. Let 

1 1 A 

y be the average of these signals and let X denote the proportion estimate associated with y. 

A 

Then X is taken as the estimate of proportions of materials in the region represented ^by aggre- 
gate of pixels corresponding to the signals y^, . . . , y^ . The components of (N^/N) X are 

taken as estimates of the proportions of the materials in the region covered by all N pixels. One 

of the advantages of averaging is that only a single proportion estimate is required for the re- 

2 2 

gion represented by many signals. Using d instead of 6 for alien object detection in the averag' 
mg procedure would require, in effect, estimation of a proportion vector for each individual sig- 
nal and thus would cause the loss of the speed advantage inherent in averaging. 
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3.4 SIGNATURE ANALYSIS 

The quality of the estimates of proportions one can expect can be determined to a large ex- 
tent by examining the pure signature set. In recognition processing we know that the quality of 
results depends upon the distances between pairs of signature means relative to their spreads 
(covariances). When these distances are large, good results can be expected. Not only is this 
requirement necessary for good proportion estimates, but a more stringent condition must be 
satisfied: that no pure signature be close in a probability sense to any signature of a mixture 
of the other materials. 

A feature of the improved proportion estimation program is a simple test called geometric 
signature analysis. We deal with the transformed signature simplex with vertices A., 1 s i < m, 
and assume that the common covariance matrix of all the transformed signatures is the identity. 
Let r^ be the distance of A^ to the closest point in the face of the signature simplex opposite A.. 
The face opposite A. is the convex hull of all the vertices A^ except for A^, Then r. measures 
the smallest distance, in standard deviation units, of A. to the mean of a mixture of the other 
materials in the signature set. If r. is small, we would expect data points representing A. to be 
confused with data points representing mixtures of the other materials. Figure 1 is an example 
of a signature simplex well-conditioned for proportion estimation. The circles at the vertices 
indicate the spread oi the distributions at the vertices; these circles are formed by points which 
are one standard deviation away from the vertex. Each vertex is several standard units away 
from the closest point in the opposite face. Figure 2, on the other hand, is an example of an ill- 
conditioned signature simplex. The pure signature mean is less than a standard deviation 
away from the closest point in the opposite face. 
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4 

DATA AVERAGING 

Two techniques designed primarily to increase speed ol calculation have been examined 
theoretically and numerically with simulated data. These are (1) data averaging, reported in 
this section, and (2) the use of a simplified estimator, reported in Section 5. 

Data averaging can be applied to problems which require an estimate of proportions for a 
larger area — say for a section or quarter -section of land, or for the land belonging to a particu- 
lar farm. In addition, data averaging may lead, under certain conditions, to a more accurate 
estimate of proportions over a wide area than does the pixel -by -pixel estimate. However, data 
averaging has no role when we desire to know the contents of the scene on a pixel-by-pixel basis. 


4.1. ANALYSIS 


Consider any area over which average proportions of given materials are to be estimated 
(see Fig. 3). Suppose that there are N resolution elements in the area, and that for each there 
is a data (signal) vector y^ and a true proportion vector X^. We want to estimate the overall 
proportion vector X for the whole region from the data vectors y^, y^, . . . , y N - Then, assum- 
ing the pixels are equal in area, 


N 

N i=l 


Let X be the standard estimator corresponding to a signal y. Sometimes, for clarity, we will 
exhibit the functional dependence by the notation X(y). The usual way to estimate X is to find V 
for each pixel, then compute the estimator X of X by 


N 

*4lX 

1 i=l 1 


This ”point-by -point estimation" of X we call Method 1. 

An alternative method for estimating X is to obtain the average y of the signals 


N 



and then compute the estimate 

A A _ 

X = My) 

This procedure is "estimation with averaging," sometimes referred to as Method 2. For larger 
values of N, the computation of X is faster than the computation of ft by a factor of about N, 
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FIGURE 1. WELL -CONDITIONED SIGNATURE SIMPLEX 



FIGURE 2. ILL-CONDITIONED SIGNATURE SIMPLEX 



FIGURE 3. SCHEMATIC FOR DATA AVERAGING 
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since only one proportion estimate is required for Method 2 as opposed to N for Method 1. We 
can now compare errors in estimates of X for the two methods. 

The basis of comparison is the mean -square error of the estimate given by 

e 2 = eI!x - xll 2 

for Method 1 (point-by-point) estimation, and by 

e 2 = eIIa - xll 2 

for Method 2 (estimation with averaging) where E is the expected value operator and II • II 
represents the Euclidean norm. For each i, 1 g- i I N let 

b. = E(\.) - X. 

l i' l 

Then b. is the bias of the estimate Let 

Then b is the average bias of the X^ Assuming the y i are statistically independent, the follow- 
ing result is derived in Appendix B: 

N 

-Eft,)!! 2 * Ilbli 2 (2) 

hp 1 1 1 

It then follows that 

Ell? - All 2 a llbll 2 ' (3) 

Now X^ and E(X^) are both proportion vectors; therefore 

V 

ll£. - E$.)ll 2 ? 2 

and it follows from Eq. (2) that 

e ^e||«-e(S)II 2 s| + llbll 2 (4) 

By (3) and (4), the expected mean-square error of the average of the estimates goes to zero 
as N goes to infinity if and only if the average bias of goes to zero; and b may not go to 0 as 
N increases. 

On the other hand, it is also shown in Appendix B that 



where 
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r = max trace M. 
i i 

1 = i = m 

and (3 depends upon A but not on N. Note that this result does not depend upon the equal con- 
variance assumption. Thus the expected mean-square error of the proportion estimate for 
the average signal goes to zero as N goes to infinity. 

Since it is reasonable to expect that the average bias b will normally remain above a fixed 
level for typical pure signature sets, X would be a better estimator of X than £ for large N. How 
large N should be depends on the pure signature configuration. Some results for this will be 
given below. 

4.2 NUMERICAL RESULTS 

Using simulated multispectral data based on five agricultural crop signatures, numerical 

A — 

tests were conducted to compare the accuracies of X and X. The test set contained about 2000 
data points; test areas were of size N = 1, 27, 135, 270, 540. Details of the simulation are de- 
scribed in Appendix C. 

2 2 

Estimates of e ^ and € 9 were obtained, with results as a function of N presented in Table 1. 
Graphs of these results are presented in Fig. 4. Note that Methods 1 and 2 yield identical estimates for 
N = 1. Note also that for smaller value of N > 1, Method 1 gives more accurate results. How- 
ever, in Fig. 4, the two curves cross at about N=330 and Method 2 remains more accurate there- 
after. This agrees with the theory of Section 4.1. From the graph for Method 1 we can estimate 
that the norm-square of the average bias of the estimates X. is about 0.007, the apparent limiting 
value. 

4.3 DISCUSSION AND CONCLUSIONS 

The results of the analysis and of the numerical tests, taken at their face value, are in 
agreement with each other and lead to the same conclusion — namely that for N greater than 
some threshold value, estimation with averaging should give more accurate results than point- 
by-point estimation. The value of N for which this occurs is dependent on the particular sig- 
nature set and on the prior distribution of the proportion vectors. However, a word of caution 
is in order. This result is based on a highly theoretical situation in which a small, well es- 
tablished set of signatures define the statistical distribution of data points assumed to be un- 
correlated. In the world of real data we find that there are always alien materials present. 

The signal values are not necessarily correlated. Finally, if we attempt to average over a 
large area to overcome these problems, we may include too many materials in the averaged 
signal; ultimately, no more than n + 1 materials can be included in the signature set, where n 
is the number of spectral channels. 
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Therefore, the details of when averaging should be used and, if it is used, the question of 
how many pixels should be averaged, will have to be worked out by experience with real data. 


TABLE 1. MEAN -SQUARE ERROR IN ESTIMATING X. (N = number 
of pixels in area represented by A..) 


Method 

N=1 

N=27 

N=135 

N=270 

N=540 

1 tf) 

0.2465 

0.0145 

0.0083 

0.0074 

. .0.0073 

2 (X) 

0.2465 

0.0385 

0.0132 

0.0075 

0.0052 
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FIGURE 4. MEAN-SQUARE ERROR IN ESTIMATING X 
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5 

A SIMPLIFIED PROPORTION ESTIMATOR 

In another attempt to improve the speed of proportion estimation, a simplified estimator 
was studied. This estimator is described in Section 5.1 and compared with the standard esti- 
mator. The simplified estimator is then analyzed with respect to averaging properties. In 
Section 5.2. numerical comparisons of the two estimators, based on simulated data, are given 
with respect to speed of computation and accuracy. 

5.1 DESCRIPTION AND ANALYSIS 

In Section 3 the standard proportion estimator was defined as the minimum of 

G(X) = 1 1 y - \ W 2 
over X, subject to the constraints 


X 1 > 0 and i = 1, . . . , m ^ 

/ 1 2 . m., 

where X = (X , A. , . . . , X 

In this manner the estimate $ of X is obtained via quadratic programming. In the simplified 
estimation procedure the minimization problem is solved subject to constraint (5) but not con- 
straints (6). Solving the problem without constraints (6) is equivalent to finding the projection 
of y onto the linear hull of the signature mean vectors (vertices of the signature simplex). 

If this projection falls within the simplex, then constraints (6) are satisfied and the solution X 
is the same with or without constraints (6). If this projection falls outside the signature sim- 
plex, then the solution will contain negative components. A proportion is then obtained by 
setting these negative components equal to zero and normalizing the remaining components. 

In order to define precisely the simplified estimator X, the following definitions are nec- 
essary. Let S denote the set of all proportion vectors of dimension n. Let S A denote the pure 
signature simplex. Then S A is the set of vectors of the form AX where X is in S. Let L be the 
linear hull of S. Then L is the set of all vectors of dimension n with component sum equal to 
one. Clearly, L contains S. Let L^ denote the linear hull of S^. Then L^ is the set of points 
of the form A rf where r) is in L. For a point 77 in L, let i denote the vector obtained by setting 
all negative components of 77 equal to zero. Let p = p(t]) denote the sum of the positive compon- 
ents of 77 . Note that p > 1 because the sum of all the components of r\ is 1 . Now let 77 = - rj ■ 
Note that 77 is a proportion vector. 
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Let Py denote the orthogonal projection of the signal y onto L^, i.e., Py is the unique point 
in closest to y. Py has a representation 

Py = A 77 

for some 7? in L. This representation is unique for non-degenerate S^. If X Q is the true pro- 
portion vector of the pixel represented by the signal y then the estimator X q of is defined 
to be 

*o = ^ 

The computational advantage of the estimator X over £ is that X may be essentially obtained 
by a single projection. With regard to accuracy, its desirable properties are obscure. If Py 

~ A A 

falls within the signature simplex S A then X = X. But in this case, the computation of X using 
the ERIM-modified Theil and van de Panne algorithm would be essentially the same as the 
computation of X, so that X is not advantageous. When Py is not in S A , X and X may differ con- 
siderably. Figures 5(a) and 5(b) show the relationship between X and^ for two configurations. 

In Fig. 5(a), X = (3/4, 1/4, 0) and X = (1/2, 1/2, 0). In Figure 5(b), X = (1,0,0) and X = (0,1,0), 
the maximal difference. 

Some properties of X with respect to averaging procedures will now be examined. These 
properties are similar to the corresponding averaging properties of the estimator £ . Let 
1 s i 5 N, denote a signal from a pixel with true proportion vector X.. The region represented 
by the aggregate of the N pixels would then have true proportion vector X given by 

N 

*-iZ\ 

i=l 1 

Here we are assuming pixels of equal size. If we estimate X by averaging the estimates X. 
from the individual pixels, we obtain the pixel-by-pixel simplified estimator (Method 1). 



Let 

b. = E(X.) - X. 

1 x 1 1 

and let 6 be defined by 
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(a) X = (3/4, 1/4, 0) and £ = (1/2, 1/2, 0) 


/ v Py = An 



(b) X = (1, 0, 0) and X = (0, 1, 0) 


FIGURE 5. GEOMETRIC COMPARISON OF STANDARD AND SIMPLIFIED 

ESTIMATORS 
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— „2 
Then b is the average bias of the estimators and the mean-square error (assuming as 

usual that the y i are statistically independent) is given by 


~2 

e l 


= E 




E(X.)|| 2 + llbll 2 


(V 


It follows that 


Ell*- fils llbll 2 , W 

Since /V and E(Xj) are both in the simplex S of proportion vectors, 

11^ - EXjl 2 S 2 
and it follows from (7) that 

i" 2 = E ll£ - xl! 2 *1 + tl6l| 2 (9) 

By (8) and (9) the expected mean-square error of X goes to zero as N goes to infinity if and 
only if norm-square of the average bias of the estimates X^ goes to zero. 

Now let 



and let \ denote the simplified proportion estimate for the signal y and consider X as an esti- 
mate of X (Method 2). A proof of the following result is sketched in Appendix B: 


~2 
6 2 


= E\\X -X 


< m ^ T 
N 


where m is the number of pure signatures 


r = max trace M. 
i 1 

1 5ism 

and £ depends upon A but not N. Thus the expected mean-square error of the proportion esti- 
mate X for the average signal goes to zero as N goes to infinity. 


5.2 NUMERICAL RESULTS 

Numerical tests were conducted to compare the standard and simplified estimators. 

These tests with the simplified estimator corresponded to the tests performed with the stand- 
ard estimator reported in Section 4.2. Thus areas of size N=l, 27, 135, 270, 540 were used 
and the simplified estimates of X were obtained via Methods 1 and 2. The mean- square error 
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of the estimates as a function of N are given in Table 2. For N = 1, X and X coincide, and their 
entry in the table is the mean- square error for the simplified estimator for a single data point. 
Comparing this value, 0.3224, with the value 0.2465 for the corresponding entry in Table 1, we 
can see that for this data, the mean-square error of the simplified estimator is about a third 
more than the mean-square error for the standard estimator. 


Graphs of the results in Fig. 6 show that is going to zero as N increases, as predicted 
by the theory. It is not clear from the graph what the limiting value 1 ^ is, but this value of 
would correspond to the norm- square of the average bias b of the estimates X. . 


For the purpose of comparing the standard and simplified estimator, the graphs of Figs. 

4 and 6 are displayed together in Fig. 7. We see that e 2 and eventually coincide. This de- 
rives from the fact that for large N, the projection of the average signal y falls inside the sig- 
nature simplex, and this results in the equality of the standard and simplified estimates. 

When we compare e ^ and we see a surprising result: £j is actually less than for 
larger values of N. We cannot explain the phenomenon. We think it may stem from the par- 
ticular simulated data set employed, and we do not have confidence in the generality of this 


particular result. 


Nov/ let us compare computation time for the simplified and standard estimates. On the 
basis of estimating proportion for each of the 2000 simulated data points, the simplified esti- 
mation procedure required about two- thirds the time required for the standard estimation 
procedure. 


5.3 CONCLUSIONS 

The standard estimation is more accurate than the simplified estimator for a single pixel. 
These two estimators are identical for larger values of N when estimation with averaging 
(Method 2) is performed. When X is estimated by Method 1, we found the simplified estimator 
superior for larger values of N. The simplified estimator requires about two-thirds the com- 
putation time of the standard estimator. We feel that the great drawback of the simplified 
estimator is the lack of theoretical underpinnings. On balance, we believe that the simplified 
estimator should be shelved for the present and that proportion estimation development should 
be continued with the standard estimator. 


28 



FORMERLY WILLOW RUN LABORATORIES, THE UNIVERSITY OF MICHIGAN 



TABLE 2. MEAN-SQUARE ERROR IN ESTIMATING X WITH 
SIMPLIFIED ESTIMATOR. (N = number of pixels in 
area represented by \.) 


Method 

N=1 

N=27 

N=135 

N=270 

N=510 

MX) 

0.3224 

0.0128 

0.0049 

0.0035 

0.0027 

2<X) 

0.3224 

, 0.0395 

0.0137 

0.0075 

0.0052 
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MEAN -SQUARE ERROR 



N (Number of Pixels in Area Represented by X) 


FIGURE 6. MEAN -SQUARE ERROR IN ESTIMATING X WITH SIMPLIFIED 

ESTIMATOR 
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FIGURE 7. COMPARISON OF MEAN -SQUARE ERRORS IN ESTIMATING 
\ WITH STANDARD AND SIMPLIFIED ESTIMATORS 
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6 

ALIEN OBJECT THRESHOLD SETTING 

In Section 3.3 the improved procedures for dealing with alien objects or materials were 
described. Alien objects are those objects not represented in the pure signature set. The 
alien object test employed a generalization of the chi-square test used in recognition process- 
ing. In mixtures processing, a signal y is said to represent a pixel with alien material if 

lly - AX || 2 > 

2 

for every proportion vector X. The alien object threshold is X Q , and in this section we are 
concerned with determining how best to set the value of this parameter in a specific situation. 

In Section 3.1, numerical tests were performed to determine the optimal alien object threshold 
setting. 

6.1 NUMERICAL TESTS 

The following experiment was performed in order to assess the effect of the alien object 
threshold setting on estimation accuracy. Three -thousand simulated multispectral data points 
were generated based on five agricultural signatures: bare soil, corn, soybeans, cut alfalfa, 
and alfalfa. (The characteristics of the data points are defined in Appendix C.) The data 
points were then divided into three sets of 1000 points each. Then each set was processed 
using four of the signatures for the pure signature set, to obtain an estimate of the proportions 
of these four crops in the region represented by the thousand data points. The point-by-point 
standard estimate of the average proportion vector for the region represented by the 1000 data 
points was obtained. This vector had four components; thus the crop not in the signature set 
played the role of alien material. Then the norm- square of the difference between the region's 
estimated proportion vector and the true proportion vector was computed and averaged over 
the three regions represented by the three sets of 1000 data points. The result was taken as 
an estimate of the mean-square error for the particular alien material and alien object thresh- 
old setting chosen. The threshold setting was varied between the values 7 and 37. Graphs of 
the results are presented in Fig. 8. The graph corresponding to a particular signature set of 
four crops is identified by the crop treated as alien. Thus, corn alien means that the classes 
represented in the signature set and used for proportion estimation were bare soil, soybeans, 
cut alfalfa, and alfalfa. 

In Table 3 we list the optimal alien object threshold for each pure signature set — or, 
equivalently, for a particular alien material. We know that 27 is near-optimal for the two 
cases: (1) corn alien and (2) alfalfa alien, because a negligible number of data points are re- 
jected at this level. 
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The graphs show that a reasonable alien object threshold setting for all five cases would 
be about 15. We did not find any clear relationship between the distance of the alien signature 
mean to the corresponding signature set and the optimal threshold value. This may be because 
the simulated data points were generated with a flat distribution over the set of proportion 
vectors. As a result the average true proportion vector was about (0.25, 0.25, 0.25, 0.25) in 
each of the five cases. The geometrical relationship of the alien signature mean to the means 
of material in the signature set may distort the optimal threshold that would result if several 
different values of X were used. This, we believe, is the reason that two cases in Fig. 8 show 
no clear-cut finite optimal value. 

6.2 CONCLUSION 

For the limited tests conducted, we found a single alien-object threshold setting satisfactory 
in all five cases. This result is based on a very limited experiment. However, the procedure 
for finding the best threshold setting on simulated data may be used on training sets for real 
data — in fact, we do this in the next section. (When there are not enough training data avail- 
able, we still don’t have a completely satisfactory way to set the alien object threshold.) 


TABLE 3. OPTIMAL ALIEN OBJECT 
THRESHOLD 

Alien Class Optimal Threshold 

12 
27 
12 
16 
27 


Bare Soil 
Corn 
Soybeans 
Cut Alfalfa 
Alfalfa 
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7 

ESTIMATING WHEAT PROPORTIONS FROM ERTS DATA 

An exercise was conducted using the mixtures algorithm to determine wheat proportions 
in selected regions from ERTS-1 multispectral scanner (MSS) data. The data were collected 
over Fayette County, Illinois, on 11 June 1973. Ground truth information was collected as part 
of the CITARS (Crop Identification Technology Assessment for Remote Sensing) task, a joint 
investigation by the Earth Observations Division (EOD) of the Johnson Space Center, Purdue 
University’s Laboratory for Applications of Remote Sensing (LARS), and ERIM. 

Specifically, the objective of the exercise was to estimate by mixtures processing the pro- 
portion of wheat in each of 20 quarter-sections. Results were disappointing in that the rms 
error of the estimates was almost 85% of the known proportion of wheat averaged over these 
quarter- sections. By comparison, the rms error reported for estimates obtained by conven- 
tional recognition processing for the same quarter-sections was about 50% [9]. 

Tests on the signature set used in mixtures processing showed that the signatures of the 
individual materials were very close in a probability sense to the signatures of mixtures of 
the other materials. Thus, we should not have expected to obtain very good results. The ill- 
conditioned nature of the signature set derives to some extent from the limited number (4) of 
the ERTS-1 sensor's spectral bands, and their relatively large width. Also, the time of year 
contributed to the low contrast. Nevertheless, results of the exercise indicate that present pro- 
portion estimation procedures should be modified for these low-contrast situations. 

7.1 DATA FLOW AND SIGNATURE EXTRACTION 

Multispectral data flow during processing as well as preliminary data handling steps have 
been already described [9 J. Procedures for extracting signature statistics from training data 
are given in the same document. We will now discuss the quality of signatures obtained and 
the procedure for selecting the signature set used in estimating proportions. 

Signatures for seven object classes were obtained: wheat, corn, soybeans, clover, trees, 
pasture, and water. The number of training data points used for estimating signature param- 
eters for each of these crop types is shown in Table 4. Note that only 9 data points were 
available for estimating the pasture signature and that 17 points were used in estimating the 
clover signature. The last column of the table lists the determinant of the estimated covariance 
matrix for each object class. [The square root of the determinant is proportional to the volume 
of the unit contour ellipse of the class signature and is a measure of the spread of that signature. 
Considerable variation among the estimated signature covariance matrices may thus be expected.] 

A training set consisting of 20 quarter -sections was selected in the following fashion: All 
of the quarter -sections are in a land segment 5 miles wide by 20 miles long. The 100 sections 
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comprising this segment were numbered starting with the top left section (No. 1) and proceed- 
ing to the top right (No. 5); the next row starts at 6 on the left and goes to 10 on the right, and 
so on. With each of the 20 quarter-sections carrying the same number as the full section to 
which it belonged, the 10 odd-numbered quarter-sections in this ordering were then designated 
as belonging to the training set. 

Since the mixtures algorithm permits, at most, five signatures in the signature set when 
four channels of spectral data are available, two signatures were eliminated by the following 
method. The distance r of each object class from wheat was computed by geometric signature 
analysis and the amount, a, of each of the classes in the 10 training quarter-sections was de- 
termined. Then the two object types for which 

<*M-r/2) 

were least were excluded. Here tf>(x) is the univariate Gaussian density function. On this 
basis the five signatures retained in the signature set represented wheat, corn, soybeans, 
clover, and pasture. 

7.2 SIGNATURE SET ANALYSIS 

Analysis of the signature set retained for mixtures processing showed that the set was 
ill-conditioned — that is, some or all of the pure signature means were close, in a probability 
sense, to the mean of the signature of some mixture of the other materials. In fact, for this 
particular case, all of the pure signature means were close to the mean of the signature of 
some mixture of the other four materials. Table 5 quantifies this separation in standard devi- 
ation units for each of the materials in the signature set. Our past experience with the mix- 
tures algorithm indicates that, for satisfactory results, all the distances determined by signa- 
ture set analysis should be about 3 or more. 

7.3 SETTING THE ALIEN OBJECT THRESHOLD 

The procedure for setting the alien object threshold was to take nine specified values and 
then choose that value found to be best for the training data. The nine values initially taken 
were all too large, so much lower values were then tried; the best alien object threshold value 
for the training data was found to be 2. This value was used for processing all 20 training 
quarter-sections. 

Only five of the ten training quarter- sections were employed for setting the alien object 
threshold — namely, those numbered 25, 32, 44, 83, 86. The training quarter -sections num- 
bered 7, 19, 72, 73 were excluded because they contained less than 10% wheat, which we felt 
would bias the threshold setting to the low side. Training quarter -section No. 56 was excluded 
because the ground truth at hand was of questionable validity. 
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TABLE 4. NUMBER OF DATA POINTS USED IN ESTIMATING 
SIGNATURES AND DETERMINANTS OF ESTIMATED 
COVARIANCE MATRICES 


Object Class 

Number of 
Data Points 

Determinant of Estimated 
Signature Covariance Matri 

Wheat 

47 

265.5 

Corn 

88 

12934.3 

Soybeans 

116 

750.4 

Clover 

17 

47.4 

Trees 

180 

3.7 

Pasture 

9 

3.6 

Water 

120 

6.5 


TABLE 5. RESULTS OF SIGNATURE SET ANALYSIS. 
(Distance is measured from mean of the material 
signature to the mean of the mixture signature 
which is closest to it in a probability sense.) 

Distance 

Material (Standard Deviation Units) 


Wheat 0.22 

Corn 0.31 

Soybeans 0.24 

Clover 0.41 

Pasture 0.17 
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7.4 MIXTURES PROCESSING RESULTS 

Table 6 gives the results obtained in using mixtures processing to estimate the proportion 
of wheat for each of the 20 quarter-sections. For comparative purposes, results of conventional 
ERIM recognition processing [9] are included. With mixtures processing, the rms error for 
the estimated wheat percent was about 0.85 of the average wheat percent over the 20 quarter- 
sections — not good. Since the corresponding rms error was about 0.6 for recognition process- 
ing, it is clear that recognition processing outperformed mixtures processing. 

7.5 CONCLUSIONS 

The fact that mixtures processing gave poorer results than recognition processing in this 
exercise indicates that mixtures processing techniques require further development for low- 
contrast situations. 

The ill-conditioned nature of the signature set, as revealed by signature set analysis, indi- 
cates that the spectral information gathered during certain times of the year by the four ERTS-1 
channels is insufficient to disperse the signatures of materials enough to permit successful 
mixtures processing without modifying present procedures. 
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TABLE 6. RESULTS OF ESTIMATING WHEAT PROPORTIONS IN 20 QUARTER- 
SECTIONS BY MIXTURES PROCESSING AND CONVENTIONAL 
RECOGNITION PROCESSING 

Quarter -Section Ground Truth Estimated % Wheat Estimated % Wheat 

(ERIM No.) Wheat % Mixtures Processing Recognition Processing 


7 

0 

20.3 

13.0 

18 

8.3 

11.1 

14.0 

19 

4.4 

4.1 

2.1 

21 

28.5 

12.5 

11.8 

25 

21.9 

19.8 

12.4 

28 

22.9 

16.2 

21.7 

32 

15.5 

8.0 

8.5 

40 

2.8 

11.9 

3.0 

44 

18.4 

21.3 

18.8 

46 

0 

14.2 

17.3 

56 

26.6 

10.0 

10.7 

64 

20.5 

4.6 

10.9 

72 

9.3 

12.2 

14.9 

73 

0 

11.5 

15.2 

75 

0 

9.1 

8.0 

82 

22.4 

15.7 

14.0 

83 

17.6 

10.3 

15.3 

84 

22.5 

25.0 

35.2 

86 

12.5 

23.0 

24.2 

100 

flooded 

(0) 

18.0 

— 

Average 

12.9 

13.8 

13.9 

RMS Error 


11.0 

7.9 
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8 

CONCLUSIONS AND RECOMMENDATIONS 

Hie multispectral data processing method for estimating proportions of objects and ma- 
terials appearing within the instantaneous field of view of a sensor system has been improved 
— especially with regard to the computational speed of the basic algorithm and the associated 
alien object detection scheme. 

Although the simplified proportion estimator did not prove adequate, results obtained in 
tests of these procedures on more or less artificial data sets indicate their potential useful- 
ness. Averaging schemes to further increase processing speeds appear to be advantageous 
under certain circumstances. 

Success with applications of proportion estimation techniques on ERTS data have been re- 
ported [4, 5]. However, during certain times of the year a low-contrast situation may occur. 
In such situations, we feel some modification of the techniques will improve performance. One 
such modification would be based on the assumption that any pixel contains a very limited 
number of materials, say not more than two or three. In addition, the estimated proportion of 
a material in a pixel should perhaps be required to have some minimum value before it is 
counted. 

We believe that advances in proportion estimation now require experimentation with real 
data sets accompanied by adequate associated ground truth information. 
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Appendix A 

DETAILS OF COMPUTATIONAL IMPROVEMENTS 

The present adaptation of the Method of Theil and van de Panne is faster than the original 
version by a factor of 14. This progress is attributable to (1) a faster machine, (2) increased 
computer memory that obviates recomputation of certain quantities, and (3) discovery of pro- 
jection theorems that permit computational shortcuts. These newer theorems — which have 
enabled a number of improvements to be made in the estimation algorithm, the alien object 
test, and geometric signature analysis — are detailed below. 


A.l THEOREMS LEADING TO IMPROVEMENT 

The Method of Theil and van de Panne computes the closest point y = A \(y) to y in the 
signature simplex by computing certain projections of y onto the simplex and its subsimplices. 
Though it does not need to compute every such projection, the method starts by projecting y 
onto the hyperplane determined by the signature mean vectors A^, . . . , A m ; it then projects 
y onto hyperplanes of succeedingly lower dimensions. The projection y* of y onto the hyper- 
plane determined by A^, . . . , A m is computed by obtaining \ = (X^, . . . , \ m )T from the equa- 
tion 


} 

• HA 

HA 

• r im f 


"hi 


V 

r ml-' 

1 

, . r l 

mm 

, . 1 0 


X m, 
A ' 


id 


(A-l) 


where F.. = A*A^ and g^ = A^y. Then y* = A L \. To project y onto a lower-dimensional hyper- 
plane determined by some but not all of the A., delete the rows and columns corresponding to 
the A. not included, and solve the reduced system of equations for the A^. If the projection 
hyperplane does not include A., A^ must be omitted from the equation and set to zero. 


All these projections would ordinarily require a great deal of computation, but the theorems 
given below afford some savings in effort. 

Theorem 1. The square of the distance from y to y* is y y - A g - A. 

2 T 
Proof. D , the squared distance, = (y - y*) (y - y') 


T T T T T T t 

= (y - A \) (y - A a A) = y y - 2A 1 Ay + A 1 AA a. 


Let J be a column vector of m l's. Equation (A.l) implies 
T\ + AJ = g 


fin this appendix, denotes the i-th component of the vector A and the Aj^ are row vectors 
forming rows of the matrix A. 
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In other words, 

AA T \ = Ay - AJ 

2 2 

Substituting in the third term of D , D becomes 

rpfTirp TT TTT 

y y - 2X i Ay + X*(Ay - AJ) = y y - X Ay - AX J = y y-X g-A 

because X^J = 1. Q.E.D. 


This theorem provides a convenient way of computing the ”out-of-plane distance” in the 
alien object test described. 

Let X’ be the X vector that would be obtained by deleting a vertex A fe from the set of vertices 
considered. The next theorem shows a relationship between X and X'. 

and C = T~ l 

Because T + is symmetrical, C is symmetrical. The solution to Eq. (A-l) is 


Let = 


r j 

J T 0 



ik 

Theorem 2. x! = X. - 7 =r-X v 

1 1 u i,i, K 


A' = A 



(A-2) 


where the subscript + stands for "m + 1”. 

Proof. Let C' be the result of inverting a that is missing the k-th row and column. It is a 
matrix theorem that 


C.. C.. 
ik jk 


C.. r = C..- n 
‘J 1] c kk 

This lemma can be verified by multiplying row i of C' by column h of r + \ The result is 

c.^X _ c„ 

: kk 


L V *V ■ £ C « - N* '2^ ' ^£ Clfcr+ih 

j^k j \ / j j 


J 


J 


The term for j = k can be included because the factor in parentheses is 0 when j - k. 
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£C.. r, ., =0 because C., = G for k # h. and C is the inverse of I\. Thus the second term 
jk +jh jk k] + 

di'ops out. The first term is 1 or 0 for h equal or not equal to i, proving the lemma. 

Define g m+1 to be 1. By Eq. (A-2) 

m+1 

x.’ = V s c.. T g. 

i L-i v J 

M 

with the k-th term missing. Therefore 
m+1 / C.,C,\ 

s' ■ z: h - -e k 


with the k-th term included because, as was seen, it is 0. This equation implies 
m+1 C.. m+1 C., 


iini . . i n ' x v . a 

QED - 


j=l KK . =1 KK 

The proof of the theorem for A is analogous. 

Because C and C’ kk depend on the A. values but not on y, they are constants that can be 
precomputed and stored; Theorem 2 provides a quick way to proceed from X to \\ thus making 
possible a considerable speedup in the Thiel and van de Panne method. 

o 

Theorem 3. The squared distance d from y* to the k-th face of the simplex is 


By the k-th face of the simplex is meant the plane through all the vertices A. except A 

1 K 

Proof. By Theorem 1, 


(T =y*y* - \' V - A’ 

T T T T 

y* y* = x AA X because y* = A X. Let us define X^ = 0 to keep X' a convenient dimension. ' 
g' = A'y* where A 1 is A with the k-th row missing. But let us define = A^y* to make g' the 
same dimension as X' and give it the convenient definition Ay*. The value of X ,T g’ is unchanged 
by this maneuvering, 

g' = Ay* = AA T X 
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Thus 

d 2 = \ t aa t x - x ,t aa t \ - A' = (X - X') T AA T X - A 
Now AA T X = g - AJ by Eq. (A-l). By Theorem 2, 



and 


A' 



SO 

d2 = c^ (cik)T<g ' AJ) + c+k ' 4 = + c+k ' 4 ? Clk ] ■ 4 

Y] C.. = 0 because it is the k-th row of C times the last column of r + . £ C^g. + C +k = X k , 

I ^ i 

by Eq. (A-2), hence 
2 \ 2 

a = -pr A Q.E.D. 

c kk 

If X. is negative, y* is on the outside of the simplex, i.e., the k-th face separates y* from 

K 

the simplex. If every X R is non-negative, y* is in the simplex. We therefore define a lower 
bound 6 of the distance from y* to the simplex as follows. 


2 /V N 

6^ = min 

\ <0 \ C kk 


- A 


if some X^ < 0=0 


if all X k 's g 0 


6 is a lower bound because the nearest point on the simplex to a point outside it has to be on some 
face, say the k-th face of the simplex. The distance from y* to that point g the distance from y* 
to that face = S. 

Corollary 4. A lower bound for the square of the distance from y to the simplex is 


Proof. Let s be any point on the simplex. Because y-y* is orthogonal to y* -s, the squared 

T T 

distance from y to s is the squared distance from y to y*, which by Theorem lisy y-X g-A, 

plus the squared distance from y* to s, which has been shown to have the lower bound 6 . 

Q.E.D. 
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The next theorem shows conditions under which the A term in the distance expressions 
can be ignored. 

Theorem 5. When either y or the origin of the space is in the plane through the vertices, 

A = 0. 

Proof . If y is in the plane of the simplex, y = Y'. X^A^. 

k 

Si = ^kV T = IVik 

m m 

a = c + T c. g. = c + 71 c. Y\ x, r., 

++ i+ & i ++ b-i i+, k lk 

i=l i=l k 

Reversing the order of summation, 

m m 

‘■'-•SAgv* 

m+1 

C is symmetric, and, Y\ C. T.. is the inner product of the m+1 row of C with the k-th row 
i=l 1+ 1K 

of r + . Hence when the last term of the sum is deleted 
m 

L c i+ r lk - - c ++ 

m 

s =c^Ev-g=o 

because the \'s sum to 1, and the first part of the theorem is proved. 

To prove the second part of the theorem, we draw a distinction between the ’’mixture space," 
which is all points of the form 

11 mm 

where + . . . + X^ = 1, though not necessarily all positive, and the ’’vector space" spanned 
by the vectors A^, . . . , A^. The mixture space is what we have previously called "the plane 
through the vertices." 

We will first show that when the origin is in the mixture space, the mixture space and the 
vector space are identical. Every mixture is a linear combination of the basis vectors and 
hence is in the vector space. It remains to show that every point in the vector space is a mix- 
ture. Let the origin be the mixture X^A^ + . . . + X m A m and let V = CjAj + . . , + c m A m be 
an arbitrary point in the vector space. Let dbeCj+... + C m . 
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V = CjA, + . . . + C m A m + (1 - dXXjAj + . . . + X m A m ) 

because the second term is 0. 

v = [Cj + (1 - dJXjJAj + . . . + [C m + (1 - d)X m ]A m 

The coordinates now add up to 1 and V is therefore a mixture. 

Now let y be an arbitrary point in the whole space, y = y* + y' where y* is in the mixture 
space and y' is orthogonal to it. 

y = (S \V T + yr 

g. = A.y = A.SX k Aj + A.y’ = EA^ 

The second term drops out because y' is orthogonal to every vector in the mixture space. The 
remainder of the proof is the same as for the first part of the theorem. Q.E.D. 

Corollary 6. The distance from a vertex to the opposite face of the simplex is l/V^kk' 

Proof . By Theorem T, the squared distance from A k to the k-th face is 



\ k is 1, and by Theorem 5, A = 0. Hence the corollary. Q.E.D. 

The application and usefulness of these theorems are illustrated in the remainder of this 
appendix. 

A. 2 DESCRIPTION OF THE ALGORITHM 

The quadratic programming problem derived from maximum likelihood and the mixtures 
model is solved by an adaptation of the Method of Theil and van de Panne. Since our previous 
report, certain improvements have been made in the algorithm. For the sake of completeness, 
our original version [ 2] of Theil and van de Panne is outlined here. A description is then given 
for the improvements which have permitted reduction of computation time. 

Let {h} denote the set consisting of the integer h, and let S be a subset of the index set 
{1, a, . , . , m}. Let the A, be the vertices of the (transformed) signature simplex. Let H S 

J s s 

be the hyperplane formed by those B. so that j < S. Let z be the projection of z onto H . 

S ^ 

Then there exists a vector p so that 



46 



Terjm 

FORMERLY WILLOW RUN LABORATORIES. THE UNIVERSITY OF MICHIGAN 


L S 

J 


= 0 


for j c S 


I 


J 


= 1 


S c 

However, X may not be a proportion vector, since there may be j so that X? < 0. The set of 
such j is called the violation set V(S) of S. Let 0 denote the empty set. Obviously, if V(S) = 0 , 

g 

then X is a proportion vector. 


As noted before, the Theil and van de Panne method finds the closest point in the simplex 

s 

to z by computing projections of z onto various hyperplanes H , in such a manner that not all 
such projections A have to be made. Let be the desired proportion vector. Then there exists 

A A S 

S such that X = X , and the following rules (theorems) apply: 

(1) If V{0 ) * 0 , there is some h e V(0 ) such that h e S. 

(2) If V(S) *0 , then SCS only if there is at least one h e V(S) so that h c S. 

g 

(3) For any S such that V(S) = 0, X is the desired proportion vector only if h e V(S - {h}) 
for all h € S 


Given these rules, it is possible to define the algorithm in detail. 


The Theil and van de Panne method is divided into, at most, m stages. In stage k, the sets 
S^. all have exactly k elements. The algorithm proceeds as follows: 

Stage 0: Compute X^. If V(<f>) = <f>, then X^ is the desired proportion vector and the algo- 

A 

nthm terminates. Otherwise, S * 0 . Proceed to Stage 1. 

S 1 

Stage 1: Using rule (1) above, consider all S.^ = {h} so that h e V(0). Compute X and 

Si 

VfSj) for all such Sj. If for any S^, V(Sj) = 0 , apply rule (3) to determine whether X is the 
desired solution. If it is, the algorithm terminates. If not, Sj is removed from the collection 
of sets for Stage 1. 


Stage k: Using rule (2), construct from the remaining sets S^, the sets for phase k 
by adding to each S k l an element of V(S k _ 1 ) to obtain one or more S k from each j. Elimi- 

^k 

nate duplicate S sets as they occur. As in phase 1, compute X and V(S , ) for each S. . Using 

S k k S. 

rule (3), stop if a solution X is found. As before, remove any set S k sothat V(S k ) = <f> and X is not 

the solution. If necessary, proceed to phase k+1. 


In the improved Theil and van de Panne the stages are renumbered, so that "stage" cor- 
responds to the number of vertices determining the projection hyperplane — rather than the 
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number of them omitted from this determination. The "unconstrained projection" is deter- 
mined at Stage m, which is still done first. 

One major improvement takes advantage of the larger memory of the IBM 7094: Inverses 
of T + and its required submatrices are all computed and stored before the start of the print- 
by-print estimation of proportions. Thus, a matrix inverse need no longer be computed for 
every projection. The \'s at Stage k - 1 can be computed from one of. these inverses and the 
corresponding \ at Stage k, by Theorem 2. These two considerations make for savings in 
computation time. 

The theorems presented above also provide for some simplification in geometric signature 
analysis (Theorem 6) and the alien object test (Theorems 1, 3, and 6). 

The present version of Program MIXMAP was written for the IBM 7094. An explanation 
of this program is given below. This includes the present computational procedures for geo- 
metric signature analysis and the alien object test. 

A. 3 A DETAILED EXPLANATION OF THE PROGRAM 

In Step 1, the Boolean variable FIRST is set = IB before the statement to read the control 
data. In Step 2, a long conditional branch starting at "WHENEVER FIRST" reads signatures, 
makes transformations, and calculates inverses of submatrices of GAMMA. In this loop, 

FIRST is set equal to OB. Thus, the branch is entered after the first PROCESS input following 
control input and at no other time. 

The entry TRUBL is usually called in POINT modules in case of difficulty. It calls 
SYSTEM if ITEST = 0 and otherwise calls ERROR, producing an octal dump. In MIXMAP, the 
function name variable TRUB is set equal to TRUBL., and then TRUB is called in various con- 
tingencies. The reason for this roundabout procedure is that if an unexplained bug appears in 
the program, you can include after the statement "TRUB=TRUBL." a statement "WHENEVER 
ITEST. G.l, ERROR2." and get a decimal dump to help find what's wrong. A decimal dump 
could, in that case, be forced by giving NSIG too large a value in the control input. A binary 
deck for XDUMP is also needed. 

The transformations made in the "WHENEVER FIRST" branch of Step 2 are as follows: 

The vertices AA(NSIG by ND) and the covariance matrices MM(NSIG by ND by ND) are 
read by READMX. ND is the number of channels on the data tape, presumed to agree with the 
number of channels in the signature decks. 

The L5 loop forms a new vertex matrix A and cumulates ABAR, the negative of the average 
vertex, and MBAR, the average covariance matrix, using only the subset of channels chosen. 
ABAR is then added to each vertex because the new origin of the space is going to be the average 
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vertex. There are two reasons for making this translation. One is to reduce A to 0, as shown 
in Theorem 6, so that A does not have to be calculated and can be left out of the in-plane and 
out-of-plane distance formulas. The other reason is to ensure an easily invertible GAMMA 
matrix. 

MXSCL. ( FSIG, ABAR, A BAR, NX) multiplies the vector ABAR (1) . . . ABAR(NX) by the 
scalar FSIG, = l./NSIG, and stores the results in ABAR. 

A covariance-removing transformation y = Px is calculated as follows. Let M be the 

T 

covariance matrix of x. It is a theorem that the covariance matrix of y is PMP . 

Proof, Cov y = £yy^ - CySy^ 
where stands for ’’expectation of.” 

Syy = CyEy = SPxx P - SPxSx A P A 
P is a constant and comes outside the expectation operator. So 

covy = P(Exx T )P T - PExSx T P T = P(Cxx T - ExEx T )P T = P M P T Q.E.D. 

The problem is to find P such that 
PMP T = I 

-1 T-l 
M = P P 

-1 T 
M = P A P 

because the inverse of a product is the product of the inverses in reverse order. Thus P is 
found by inverting M, making the Cliolesky decomposition 

-1 T 
M = P P 

and then transposing P. 

In the program, MBAR is taken as the best estimate of M, then inverted by GJR with the 
result stored in MBAR, and then subjected to the Cholesky decomposition subroutine LLT with 
the result stored in MBAR. Because the result of LLT is a lower triangular matrix, all ele- 
ments above the diagonal are zeroed and removed. The transpose of MBAR, an upper tri- 
angular matrix, is the transformation to apply to the data vector X. 

However, to save time in transforming X, the multiplications by zero are omitted. In the 
double loop early in PROMIX, the point-processing internal function, the last value of y, 
namely Y(NX) is calculated first. X(NX) is the floating equivalent of the input integer 
DATUM(NX) and Y is the initial value of Y plusX(NX) times the (NX, NX) value of the transposed 
MBAR, now called P. Y(NX-l) depends on X{NX), X(NX-l) and two values of P. Doing the loop 
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backwards this way allows the X’s to be calculated as they are used. The initial values of Y 
are ABAR, the negative of the average vertex, thus effecting a translation to put the origin at the 
average vertex. The P's are referenced by the linear subscript K to avoid the time to compute a 
double subscript. Thus P is stored, not in matrix form, but in the unusual linear form required 
by the backwards loop in PROMIX. 

The vertices A and the average vertex ABAR must also undergo the same covariance - 

removing transformation. Any column vector Z would be transformed PZ, but if Z were a row 

vector, it would be transformed ZP . The vertices are stored in A as row vectors and ABAR 

can be considered a row vector. Thus A and ABAR are transformed by multiplying them by 

T 

MBAR, which has been reduced to P . 

The GAMMA matrix is computed next. GAMMA(I,J) is the inner product of the I-th and 
J-th rows of A. A column of l's and a row of l's are added, and GAMMA{NSIG+1, NSIG+1) is 
defined as 0. We have found the GAMMA matrix often poorly conditioned, that is, attempted 
inversion of GAMMA led to trouble. The problem was solved by making, in effect, a change 
of scale 

new Y = old Y/V 

where V is the square root of the largest element W of GAMMA in absolute value. The effect 
on GAMMA is to divide the first NS1G rows and columns by W, with the result that the largest 
element is now 1. A, ABAR, and P are also divided by V, so that the transformation takes no 
time at the point-processing stage. 

The L20 loop goes through every subset of the vertices, computing and storing the inverse 
of the corresponding submatrix of GAMMA. A subset is identified by an integer ISUB which, 
when looked at in binary form, has a 1 corresponding to every retained vertex. A subset of 
the first, second, and fifth vertex, from right to left, has the number 10011 in binary, which is 
the same as 19 in decimal. Every subset of five vertices can be represented as an integer 
from 0 to 31. The L20 loop eliminates subsets of size 0 and 1 because the corresponding in- 
verses are not needed for the projections. 

Within the L20 loop, the subset is expanded into the vector ICOL, so that ICOL(I) = 1 if the 
vertex is present and 0 otherwise. The same vector is also called BCOL so that it can be 
treated as a Boolean variable as well as an integer. The sum of the vector is the number of 
vertices in the subset, and is stored in ICOL(O). This number plus 1 is called ICOL1 and is 
the row length of the corresponding submatrix of GAMMA, It is one more because of the 
border of l’s. ICOL{NSIG+l) is kept equal to 1 to indicate that this border is always part of 
the subset. The rows and columns of GAMMA are picked out where ICOL(I) = 1 and stored in 
C(H-+T) . . . C(H+ICOLl*ICOLl). The inverse is then computed using subroutine GJR and 
stored in the same place. The index H starts at 0 and after every trip through the loop it is 
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incremented by ICOL*ICOLl. The reason it is not incremented by ICOL1 2 is that the last row 
of the inverse is not used in calculating the projections and therefore does not have to be 
saved. 

To retrieve the inverses that have been computed and stored sequentially in the C vector, 
a pointer CLOC(ISUB) is saved. Logically, the starting index H should be saved by CLOC. But 
^ - ICOL1 is saved instead, so that the zero-th element of the k-th row of the inverse is the 
simple expression 

C(CLOC(ISUB) + K*ICOLl) 

The distance from each vertex to the opposite face can now be obtained from the inverse 

C of the subset NSUB, the one that includes all the vertices. By Theorem 4, this distance 

D(K) is for each vertex K. To restore it to the original standard deviation units, it is 

multiplied by V and printed out as DIST(K). The index of is 

KK 

CLOC(NSUB) + K*(NSIG+1) + K 
which reduces to 

CLOC(NSUB) + K*(NSIG+2) 

DMIN, the smallest of these distances, is a measure of how poor the estimation of proportions 
is likely to be. It has been stated that if DMIN were 0, the estimates would be totally ambiguous. 
If DMIN for a certain K were near 0, then small chance fluctuations in X(K) would cause wild 1 
fluctuations in the mixture estimates. The program limit, MINDIS, has a default value of 0.5. 
Many users, however, may want a more stringent limit, and require, for example, that all 
vertices be at least one standard deviation away from the opposite face. Then they would set 
MINDIS = 1 in the control input. Any DMIN less than MINDIS causes the program to stop. 

Some Boolean switches are precalculated to speed up the point-processing. BAL means 
"make the alien object test," BOUT means "compute the out-of-plane distance," BAV means 
"keep a cumulative sum of the data points" so that mixture estimates can be obtained for the 
average data point, NOMIX means "do not compute mixture estimates point by point," NOXAL 
means "compute neither mixture estimates nor alien object tests," and ITEST2 and ITEST3 
signal a debugging printout. The average point AVEPT (an integer variable), PRO (the cumu- 
lative proportion vector), ALNO (the number of alien objects), and HNO (the number of non- 
alien objects) are zeroed prior to point-processing. 

The point -processing routine is the internal function PROMIX. "WHENEVER NOXAL, 
TRANSFER TO UPDATE" is used because if there are neither mixture estimates nor alien 
object tests to calculate, all that remains is to update AVEPT and HNO and return. 
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The procedure for getting Y has already been described. The floating point equivalent of 
the integer data vector DATUM is stored in X. Y is computed with the 1/V transformation 
built into P, and a translation of coordinates to move the origin to the average vertex is ef- 
fected. The vector G is calculated. G(K) is the inner product of the K-th vertex with Y. 

The initial A values are calculated — that is, the A's of the projection onto the vertex 
space. A is just the inner product of the K-th row of CL, the inverse of GAMMA, with G. If 
all the A’s are positive, the estimation is complete, and we can transfer to OUT. ICOL, L, 

SETL and STM are given values first, because these variables are used in the section of pro- 
gram starting with OUT. 

In the L110 loop, MINL is computed to be min |LAM{K) *D(K)} or 0, whichever is less. 
D(K), we remember, is 1 Thus by Theorems 3 and 6, MINL^ is the square of the in- 

plane distance for the alien object test. The out-of-plane distance squared is, by Theorems 1 
and 6, the inner product of Y with itself minus the inner product of LAM with G. D2 is the 
sum of the two squared distances unless no out-of-plane distance is calculated, in which case 
it is only the square of the in-plane distance. 

Theoretically, one would take the square root of D2 and compare it with ALSIG, which is 
the distance in standard deviations Y has to be from the simplex to be considered alien. ^D2 
is actually a lower bound for the distance. The out-of-plane distance is exact, but the in-plane 
distance is the largest distance from the point to a face that separates the simplex from the 
point, and is therefore a lower bound for the true in-plane distance. Thus, if A (D2 > ALSIG, 
then the true distance > ALSIG, the point is alien, and the mixture estimates need not be com- 
puted. If Vd 3 < ALSIG, then proportions can be calculated, a true distance obtained, and the 
test reapplied. Thus the alien points, and only they, are rejected.' In practice, a precomputed 
value TEST=ALSIG /W is compared with D2 to allow for the 1/V transformation and to avoid 
taking the square root. 

When a point fails the alien object test, the alien object counter is updated, DATUM(NSIG+2) 
is set = 1 (indicating the point is alien), DATUM(NSIG+1) is set = D2, and DATUM(l) . . . 
DATUM(NSIG) are given the conventional value 777g. 

When a point passes the alien object test and the option of only calculating mixture esti- 
mates from the average point has been selected, then there is nothing left to do but transfer to 
UPDATE, where the average point and the good point counter are updated and the function re- 
turns. 

The two-page TVLOOP carries out the Thiel and van de Panne procedure with a few modi- 
fications. The loop passes in stages from NSIG to 2. At each stage ST, it is assumed that the 
A’s, namely LAM(J), have been computed and stored sequentially, and the action at stage ST is 
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to compute the X’s for the next lower stage, and while doing that to test signs to see if all X s 
are positive. If so, the optimality test is performed, and if this test succeeds, then the esti- 
mation procedure is over. 

X values from only two stages need to be kept in storage: those at stage ST that are being 
referenced, and those at the next lower stage that are being computed. The LAM array where 
these X’s are stored is in two halves. The X's computed for stage NSIG before entering the 
TV LOOP are stored starting at LAM(O) while those computed at stage NSIG for stage NSIG-1 
are stored starting at LAM(350). After stage NSIG has been completed, and all X's computed 
for stage NSIG-1 (and have been tested and found wanting) then the stage NSIG X's are no longer 
needed and at stage NSIG-1, the X’s for stage NSIG-2 are computed and stored starting at 
LAM(O). 

At stage ST the sets of X’s that have been computed are gone through using the index J. 

J starts at Jl, which is initially 0 and is subtracted from 350 after every stage. Thus J1 alter- 
nates between 0 and 350. Preceding each set of X's in the LAM vector is a word giving the 
subset to which the X’s refer. This word is given the integer name SET(J) and made equivalent 
to the floating point array LAM. Thus, at the outset SET(O) = NSUB and the initial X's are stored 
in LAM(l) . . . LAM(NSIG). 

The JLOOP going through the sets of X’s is almost as big as the TVLOOP. The increment 
is by ST+1 to allow for the defining subset SET(J) as well as the ST X’s stored compactly, with- 
out gaps. The subset is expanded into ICOL(l) . . . ICOL(NSIG), where ICOL(N) = 1 if the N-th 
vertex is present in the subset, and 0 otherwise. 

Within the JLOOP (and almost as big) is the KLOOP, going through the X^s for K = 

1 . . . ST. Actually the loop goes for N = 1 ... NSIG and only executes when ICOL(N) = 1, at 
which point K is incremented. In this way, both the index K of the X’s, which are stored com- 
pactly, and the index N of the G vector, which is not, are available for use in the KLOOP. 

We are looking within the KLOOP at a value of N for which 1C0L(N) = 1 (i.e., BCOL(N) = 

IB). SETL is the original subset SETJ with column N zeroed. This is accomplished by bit- 
wise intersection of SETJ with COL(N), where COL is an integer vector 1, 2, 4, 8, 16, ... , 
which is in binary form, 00001, 00010, 00100, 01000, etc. The sign of X R is tested. If it is 
negative, then the Thiel and van de Panne procedure designates. SETL as a possible subset on 
which the nearest point to Y might be found. Therefore X's are computed for SETL, using 
Theorem 2 and the inverse of the submatrix of GAMMA located at CLOC(SETJ).* The X's 
are stored by the index L in the half of the LAM array not used by J. This is accomplished 


^Theorem 2 states that new X. = old X. - C,. * {\ V /C VV ), where it is the K-th vertex that 

, , , , , 1 1 IK IV IVIV 

has been deleted. 
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by starting L at 350- Jl outside the JLOOP. SET{L) is given the value SETL, and the new X's 
are stored in LAM(L=1) . . . LAM(L+ST-1). Then L is incremented by ST. 

If the sign of \ is not negative, then it appears that subset SETL need no longer be con- 
sidered; if any subset of SETJ contains the nearest point to Y, then that subset will be found in 
one of the SETL’s corresponding to a negative X. Thus as far as succeeding subsets of X at 
stage ST are concerned, SETL is done and need not be looked at again. Either it was a possi- 
ble subset and X’s were calculated for it, or it was a subset thrown out because of correspond- 
ing to a positive X. Thus before the test of the sign of X R , it is recorded that SETL is by- 
passed by setting D0NE(SETL)-1B. The DONE array starts out all zero. In every subsequent 
pass through the KLOOP, when DONE(SETL) is found to be IB, then that pass through the 
KLOOP is bypassed. 

The new X’s are tested for sign as they are computed. If all are positive then the opti- 
mality test is performed. This consists of going through the empty columns I of SETL (loop 
L120), forming SETI by adding column I to SETL, locating the inverse matrix C corresponding 
to SETI, and then computing Xj directly as the inner product of the NK-th row of the inverse 
with G. NK is found by counting the l's of SETI from right to left and stopping after the 1 in 
column I. It is needed because the inverse is stored compactly. Whenever X^ > 0, the test has 
failed and a transfer to FAIL is made. If loop LI 20 is completed without failure, then the last 
X’s computed define the nearest point on the simplex to Y and a transfer to OUT is made. 

When a subset SETL fails the optimality test, then the point on SETL nearest to Y has 
been found and it is not optimal. That means that neither SETL nor any of its subsets can be 
optimal, so we might as well record them all as being done. The section of code from FAIL 
to the END OF CONDITIONAL following L125 does this job. The method is to form all subsets 
of SETL consisting of 1 vertex each, and then combine them in all possible ways, recording 
DONE=lB for each combination. (The best way to understand this section of code is to work 
it through with an example subset.) 

The TVLOOP finishes with ST =2 when X’s for subsets of size 1 are calculated. According 
to theory, some subset must have had positive X’s and passed the optimality test, because 
there is always some point on the simplex that is nearest to Y. The program allows for the 
possibility that because of roundoff or a program bug, no optimal subset can be found. A mes- 
sage is then printed giving the line and point number so that the point can be examined later. 
The point is counted as alien. If ITEST ^ 2, the run ends there with a dump. (Occurrences of 
this message would be of interest to the authors.) 

To minimize the possibility that no subset will be found optimal, the tests of sign are 
made not with zero but very small numbers: Z = 0.000001 and ZM = -0.0000005. A set of X’s 
is considered all positive if they are all >ZM rather than all >0. To fail the optimality test, 

X. has to be >Z rather than >0. 
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The section of code from OUT up to STEP(5) describes the action taken after the optimal 
subset has been found. The most recently calculated set of X's, which are stored compactly, 
are spaced out and stored in the vector PROP with zeroes stored in the spots corresponding 
to the missing vertices. The true distance squared, ^~)y. - A^A.) 2 , is calculated and com- 
pared with TEST to see for sure whether the point is alien. (The alien object test represents a 
screening operation. Every point failing the test is alien, but some points passing the test 
will also be alien.) If so, the alien object counter ALNO is updated and DATUM(NSIG+2) is 
set - 1. 

If the point passes this second alien object test, the good point counter HNO is updated 
and DATUM(NSIG+2) is set = 0. If a mixture estimate will be made from the average data 

point, then the cumulative sum of data points AVEPT(l) AVEPT(NX) is updated. 

The cumulative sum of mixture estimates PRO(l) PRO(NSIG) is updated. The squared 

distance is converted to standard deviation units by multiplying by W and then scaled into the 
range 0 . . . 511 by multiplying by 5. The last multiplication is done first in Step 2 resulting 
in a factor W5. The number of output channels NX is set = NSIG+2. The point-processing 
routine PROMIX has now completed its work so it returns. 

In Step 5, the branch of the MIXMAP module called after the points of the requested rec- 
tangle have been processed, the total number of points, the number HNO of good points, and 
ALNO of alien points, and the time to process the points are printed. If mixture estimates 
have been calculated, then the sum vector PRO of the mixture estimates is divided by the 
number of good points and the result stored in the vector PROP. Internal function PRTMIX is 
called to print the estimates. 

If mixture estimates based on an average of the data points have been requested, then the 
sum vector AVEPT of the data points is divided by the number of good points, converted to an 
integer and stored in DATUM. We need now only call PROMIX to get, in the vector PROP, 
the mixture estimates based on the average data point. We must first be sure to turn off the 
alien object test, the average point cumulation, and the no-mixture-estimate bypass. The 
distance from the average point to the simplex in standard deviation units is computed by 
taking the square root of DIST*W and is then printed. 
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Appendix B 

PROOFS OF THEOREMS RELATING TO AVERAGING PROCEDURES 

The purpose of this appendix is to derive results for use in Section 4 on averaging pro- 
cedures with the standard estimator, and/or in Section 5.1 on averaging procedures with the 
simplified proportion estimator. 

We will assume the positive integers m and n are fixed with m * n + 1; here > m represents 
the number of pure signatures and n represents the number of spectral channels. S denotes the 
fundamental m-1 simplex, i.e., the set of all points in the positive octant of Euclidean m space 
E m , with components which sum to one. Equivalently, S is the convex hull of the unit vectors. 

Material i of the m materials is assumed to have a signature with mean A. and covariance 
matrix M.. A will denote the matrix which has A. for its i-th column. Let A be a column vec- 
tor representing a mixture of materials, i.e., the j-th component A^ of A is the proportion of 
material j, 1 s j < m. It is clear that A is in S. Such a vector is termed a propoition vector. 
The ERIM model for signatures of mixtures states that the mean and covariance matrices 
M. of mixture of materials with proportion vector A are given by: 

A 

A x =AX 

m . 

M x =£ M i x 

H 

Let y denote a signal obtained from a pixel containing materials with associated true proportion 
vector A q (unknown). The standard estimator A q of A q is defined as follows. The set of points 

AA A in S 

form a simplex S. which is called the signature simplex. Let Z be the point in which is 

A 

closest (in the Euclidean metric) to y. Then Z can be represented as 

z = a ! 

o 

where is a proportion vector. A^ is taken to be an estimate of the true proportion vector 

A q . is unique if the signature simplex is nondegenerate, i.e., if S A has positive (m-l)-di- 

mensional volume. We will always assume this condition is satisfied. If the pure signatures 

are Gaussian and their covariance matrices are equal, then after an appropriate preprocessing 

transformation of the pure signatures and data points, A^ is the maximum likelihood estimate 

of A . 
o 

Let us consider the following situation. Suppose there are N data points y. f lSi§N, each 
representing a pixel with true proportion vector A.. The region represented by the aggregate 
of the N pixels would then have true proportion A given by 
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X 



Here we are assuming pixels of equal size. One may estimate a proportion vector X. for 
each i-th pixel with signal y.. Then the average of these estimates X. given by 



1 


may be taken as an estimate of X. The mean square error of the estimate is then given by 

eIIx-xII 2 

where E is the expected value operator and || - j| represents the Euclidean norm. For each 
i, 1 ;S i 1 N.let 

= E(X.) - X. 

Then b. is the bias of the estimate X.. Let 



1 


Then bis the average bias of the X.. Assuming the y. are statistically independent, the follow- 
ing result may be derived. 

Theorem 1: Ei|\ - x|| 2 = -^V'eIIx - E(X.)|| 2 + ||b|| 2 

1 1 

Proof: 

eIIx -"xl 1 2 = e i lx - e(x) + E(X) - x 1 1 2 = E ! I X - E(£) 1 1 2 + 1 1 E$) - x 1 1 2 
= ^ E ||Z?i ' ^IJ 2 + ^ 2 Q-e.d. 

An alternative method for estimating X is the following. Let 



1 


denote the average of the signals. Treat y as a signal from a pixel and estimate a proportion 
vector X and take X as the estimate of X. Then the mean-square error of X given by 
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EiU - \W 2 iff 

satisfies (still under the assumption that the y\ are statistically independent), 

Theorem 2: Efix-x[| 2 s§I 

N = 

where 

t = max trace M. 

1 

1 

I g i < m . . 

and p depends upon A but not N. 

The proof of Theorem 2 depends upon three lemmas. 

Lemma 1: Let y be a data point and let X be the standard estimate of the true proportion 
vector X corresponding to y. Then 

II AX - AX II § 1 1 y - A\|| 

Let L denote the linear hull of the fundamental simplex S; i.e., L is the set of vectors with 
components summing to one. Clearly, L contains S, the fundamental simplex. The linear hull 
L a of the signature simplex is defined to be the set of points in signal space of the form 

Arj for 7] in L 

Lemma 2: Let L be the linear hull of S with A nondegenerate. Then a p (depending on A) 
such that 

111? - ^ ! 1 2 ^ p 1 1 A7 ? - All! 2 

for all 7) and | in L. (p may be taken as the reciprocal of the smallest eigenvalue of B T B 

where the (n + 1) x (m) matrix B is obtained from A by the addition of a row of l r s. The fact 
T 

that B B is nonsingular is a consequence of the fact that the signature simplex is 
nondegenerate.) 

Lemma 3: Let ML, 11 il m be positive definite symmetric matrices. Let y be a vector- 
valued random variable with mean zero and covariance matrix M • given by 

A 

m 

m * 

H 

where X is in S and X 1 denotes its j-th component. Then 

Elly!! 2 i r 
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where 

t - max trace M. 
i 1 

lli^m 

Proof of Lemma 1: Consider the triangle with vertices y, AX, and AX. Let a denote the 

" A ' A 

angle at vertex AX. Since AX and AX are both in S^, every point on the line segment between 
them is in S^. If a < ^ radians, there would be a point AX q on this segment closer to y than 
AX. This would contradict the choice of X. Therefore a ? ^ radians, and consequently the 
length of the largest side of the triangle is J J y - Axl). Since II AX - AX 1 1 is the length of 
another side of this same triangle, the lemma follows. 

Proof of Lemma 2: Let B be the matrix obtained from A by adding a column of l’s. 
Since A is nondegenerate, 

T 

B B is nonsingular 

Let r) and | be in L and let u> = t) - Then the sum of the components of tn is zero and it 
follows that 

T T 

A Aw = B Bu> 

Then 

1 1 A?? - A£|)^ = 1 1 Aco 1 1 2 = <A T Aco, oi> = <B T Bu>, co> § ck [ i col 1 2 

T 

where a is the smallest eigenvalue of B B. This last equality is a well known property of 
positive definite symmetric matrices. Thus 

ll»7 - III 2 = llcol 1 2 S^IIAt? - A^ll 2 Q.E.D. 

Proof of Lemma 3: 


m m . m 

= trace = trace = ^X J (trace M^) 5 t^Px j = r Q.E.D. 


Proof of Theorem 2: There are positive numbers J3 and r such that 
Ellx - III 2 = 1 1 AX - Axil 2 (by Lemma 2) 


= 0E I iy - AX| | 2 (by Lemma 1) 


59 




(Equation Continued) 



(by Lemma 3 ) Q.E.D. 

We now turn our attention to theorems concerning averaging procedures using the sim- 
plified proportion estimator X. As before, let L denote the linear hull of the fundamental 
simplex S and let L A denote the linear hull of the signature simplex S A - Recall that for a 
point ?j in L, 77® denotes the vector obtained by setting all negative components of V equal to 
zero. Let p = p(r/) denote the sum of the positive components of rj. Note that p 21 because 
the sum of all the components of V is 1 . Now let 77 = ^77® Note that 77 is a proportion vector. 

Let Py denote the orthogonal projection of the signal y onto L^, i.e., Py is the unique 

point in L. closest to y. Py has a representation 
A 

Py = A77 

for some 77 in L. This representation is unique for nondegenerate S A * If \ Q is the true 

proportion vector of the pixel represented by the signal y, then the simplified estimator 

X of X is defined to be 
o o 



Some properties of X with respect to averaging procedures will be examined. These 
properties are similar to the corresponding averaging properties of the estimator X. As be- 
fore, let y. , II i S N, denote a signal from a pixel with true proportion vector X^. The region 
represented by the aggregate of the N pixels would then have true proportion vector X given by 



Here we are assuming pixels of equal size. If we estimate Xby averaging the estimates X. 
from the individual pixels, we obtain the estimate 

X. = ^Vx. 

1 N/ , 1 

i=l 

Let 

b. = E(X.) - X. 

1 1 1 
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and let b be defined by 

b ^L b i 

l 

Then b is the average bias of the estimators X^, and the mean-square error of the estimator 
X is given by 

eIIx -"xll 2 

Assuming the y^ to be statistically independent, we have Theorem 3: 

Ellx - x || 2 = ^V n ||(X. - E(X.)|| 2 + llbll 2 

N z ^-i 1 1 

Proof: 

Ejjx - x il 2 = e !ix - e(x) + e(x) -xll 2 = e llx - e(x )|| 2 + IIe(x)-x || 2 

= E ||^i‘ E( V)|| 2+ flb|[ 2 =-^J]E!!x i - E(X.)|| 2 + ||b |[ 2 Q.E.D. 

Now let 

y = ^i 

and let X denote the simplified proportion estimate for the signal y and consider X as an 
estimate of X. Then 

Theorem 4: eIIx-aIIs 5 ^ 

where m is the number of pure signatures 

t ~ max trace M. 
i 

1 g i < rn 

and |3 depends upon A but not N. 

The proof of this theorem requires the following additional lemma: 

Lemma 4; For all rj in L and X in S 

II T? - xll 2 ^ mil r? - xll 2 
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Proof of Lemma 4: Let r\ e L and X 6 S. II?)eS, then ?? = € and the result follows. 
Therefore, we assume that rj t S. Therefore V has some negative components. Let u denote 
the number of negative components of ??. Then since rj also has positive components, it 
follows that 

lSi;lm-l (B-l) 


Let v° denote the vector obtained by setting the non-negative components of rj equal to 
zero. Then 


© © 

n = v + v 

0 © 

r\ is orthogonal to ij 


(B-2) 

(B-3) 


Let Xj denote the vector obtained from X by setting the j-th, 1 = j = m, component of X 
equal to zero when the corresponding j-th component of r) is negative. Let X 2 denote the vec- 
tor obtained from \ by setting the j-th, 1 ^ j < m, component of X equal to zero when the cor- 
responding j-th component of rj is non-negative. Then 

(B-4) 


\ = +X 2 


X^is orthogonal to X 2 

IU 1 II 2 a l 
UJI 2 s l 


(B— 5) 
(B-6) 
(B-7) 


Q 

is orthogonal to v 

© 

\ 2 is orthogonal to V 

Let o be the sum of the positive components of v and let p = 1 + e. Then 
€ > 0 
Now, 



(B-8) 

(B-9) 


(B-10) 


(by def. of v) 


(by Eq. B-4) 


62 



Jm 


FORMERLY WILLOW RUN LABORATORIES. THE UNIVERSITY OF MICHIGAN 




+ W\ 2 W 


(by Eqs. B-5 and B-9) 


d + 0' 


7 ] - Xj - eX^ + 1 1 A_2 


< | [?r - X a - €X 1 | | + Ilx 2 lr 
“ 0 | 1 ,e " X l|| + €llX l ll ] 2 + IIX 2 1|: 


r, @ -X. 


+ 2ellx l ll* t?® - X 1 + € Z |ix 1 11* + llXgll 


7?® - Xj + 2e - X 1 + e* + |(X 2 II 


.© 


Therefore, 


II 5 - x tl 2 § 77 ® - Xj + 


2e 


T) & - X. 


+ e 2 + llxJI 2 


(by Eq. B- 6 ) 


(B-ll) 


Now consider II 7]® || 2 . The sum of the components of 77 ® is -e. Therefore, the minimum 

1 1 (b 1 , 9 G -e 

value of II ?? II occurs when 77 has each of its v nonzero components equal to — . Thus, 


\U Q W 2 i< 


and from Eq. (B-l), it follows that 
2 


Now, 


1 

II 77 - xll 2 = 1 1 77® + v Q - \ - x 2 J| 


T} & - X. 


-X, 


(B— 12 ) 

(by Eqs. B-2 and B-4) 
(by Eqs. B-3, B-5, B- 8 , B-9) 


v - X, 


© . 

V - X. 


S Hn® - A 


+ lh e || 2 -2<>7 y , X 2 > + I lx, 

2 m n2 


\V 0! + 


+ -W + llxJI 2 


m - 1 


(by Eq. B-12) 
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Therefore, 


\\v - aip ^ 

Then, 


V ® - X, 


2 + 


ml! 7? - Xll 14 - II V - X 1 1 ^ ? (m - l)||»l® - X 1 

+ (m - l)lt\ 2 l| 2 


- 2€|j7?® - 




(B-13) 

1 2 
m - 1 

(by Eqs. B-ll and B-13) 

"i 2 


SO Q.E.D. 

R. Kauth has pointed out that m is not the best possible bound. He has outlined a proof that 


*7-Xlrs Clltf-X 


m 


where 


C 


m 


m)' 


4m 
m -t- 2 


m odd 
m even 


and that C is the best possible bound, 
m 

Proof of Theorem 4: Let P denote the orthogonal projection of signal space onto L^, the 
linear hull of the signature simplex S^. Let L denote the linear hull of all proportion vectors. 
Then 

Py. = A??, for some unique in L, 1 ^ i ^ N 
Set 



1 


Then Py = At; and there are positive -number 0 and t such that 
E lix -All 2 = E I [ 7? - All 2 


S mE I i T7 - A 


2 
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(Equation Continued) 

s m/3E I i A 7 ] - A\ 1 1 2 
= m/3E 1 1 Py - AX 1 1 2 

£ m/3E( 1 1 Py - A X 1 1 2 + 1 1 y - Py 1 1 2 ) 

= mf3E I iy - AX || 2 

The last equality follows from the fact that y - Py is orthogonal to L 
and AX. Continuing, we have 

E(|X - x [| 2 £ m^Elly - Ax || 2 


= m/3E 




g -2^ Q.E.D. 


(by Lemma 4) 


, which contains Py 


(by Lemma 3) 
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Appendix C 

DESCRIPTION OF MIXTURES DATA SIMULATION PROGRAM 

To facilitate testing of the ERIM proportion- estimation algorithm, a computer program 
was written to generate simulated multispectral data of the type that would come from a mix- 
ture of materials. This was necessary to insure that completely reliable "ground truth" was 
available for each data point, so that meaningful mean-square errors of estimation would be 
calculated over large areas. Written for the IBM 7090, this program is described in detail 
below. 

The mixtures data generated by Program SIM are based on the ERIM mixtures model, and 
may be either pure mixtures of known materials or may contain a proportion of alien material. 
Any signature set may be used to generate the data, including signature tapes written on the 
CDC 1604. The actual distribution of alien materials and the proportion vectors are governed 
by parameters set by the user. The program is written as an External Function in MAD on the 
IBM 7094. 

C.l THE MIXTURES DATA MODEL 

No matter whether the pixel contains a "pure" material or a mixture, it is assumed that the 
signal s from any given pixel is a Gaussian random variable with mean vector A and covariance 
matrix M. The ERIM mixtures model relates the statistics of the pure signatures to the statis- 
tics of the mixture according to the formulae 

m . 

A = mixture mean 

i=l 

m . 

M =^|a. 1 M. = mixture covariance 
i=l 

where X 1 is the proportion of the i-th "pure" material in the pixel, A. its signature mean, and M. 
its signature covariance matrix. This can be extended to the case where there are two classes of 
materials, each with its own proportions X 1 over the class. Let there be t materials in the first 
class and u materials in the second class. Let the pixel be completely covered by the two classes 
and let i; be the proportion of the second class in the pixel. Then the model implies that 

A ' ‘I - «>£> lA i + fJ> i+t Vt 

i=l i=l 

+ £> i+ ‘M i+t 

i=l i=l 
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The two classes might be taken to be "user" and "alien" materials, where "user" materials are 
known materials for which suitable signatures are available. Then £ is the proportion of alien 
material in the pixel. Note that £ and the X 1 are themselves random variables over a collection 
of pixels. 


For each pixel, we desire to "randomly" choose A and M such that they specify an appro- 
priate mixtures distribution. Before this can be done, it is necessary to choose random values 
of £ and the A 1 . Given the number of materials in a pixel, the A 1 are to be uniformly distributed 
The number of "user" signatures or "alien" signatures actually employed to generate the data 
is also a random variable and is taken to have a certain distribution over the set (1, 2, . . . , t). 
The distribution function F(x) of £ used in the program is 


( 0 for x < 0 


F(x) = l 


a + (1 - a - / 3 ) — 

1 



for 0 5 x - 1 


^1 for x & 1 


where oris the probability that the pixel contains "user" material only, /3 is the probability that 
it contains "alien" material only, and y is another parameter. 

For each pixel, the program must determine randomly how many "user" materials and how 
many "alien" materials are in it, and this requires some prior statistical information. It would 
suffice to specify the quantities 

u , A 
p. and p. 

which are the probabilities that the pixel contains exactly i "user" materials and j "alien" ma- 
terials. In an agricultural application, these quantities are related through the geometry of the 
field configuration and can be computed from the parameters r and r A . Both of these repre- 

sent an average ratio of field length to pixel edge, but one is used for the p . and the other for 

a 1 

the Pj . The following formulas were derived: 

Pj(t) = (1 - r) 2 

P 2 (t) = 2t - 2.5r 2 
p 3 ( t) = t 2 
p 4 (r) = 0.5 t 2 
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assuming that there can be no more than four fields in a pixel. These formulas can be used for 
either p A of p . In the existing program we have added (for "user’ 1 materials) 

A U, 

p 5 (t) - 0.25t 2 
p 6 (t) = 0.25 t 2 

since, for a very small field size, it might be possible to have as many as six fields overlap in 
one pixel. Program SIM allows a maximum of six user materials and three alien materials per 
pixel. For both alien and user materials, the p.( i) are normalized to sum to 1. Once these p. 
are determined, the program can randomly choose the numbers t and u of materials in each 
pixel and then ’’randomly" choose the proportions X 1 . 

Given the extended mixtures model 



the A. and M., and the random quantities £ and X , the mixture signature mean A and covariance 
M can be computed. A sample from the Gaussian distribution with mean A and covax-iance M is 
taken as the signal vector for the pixel. 

C.2 A DESCRIPTION OF EXTERNAL FUNCTION SIM 

All of the programming required to generate the data is done in SIM, which is written as an 
External Function. Writing the data out on tape in standard ERIMmultispectral format is accom- 
plished by OUT 94, a special output routine. Both are called by a small enabling program 
which first reads certain data (see below) and instructs the operator to mount a signature tape 
from which the simulated mixtures data are to be generated. The first call to SIM (from the 
enabling program) computes certain variables and sets up variables required for linkage to 
OUT94. Succeeding calls are made by OUT94 to Internal Function LINE (part of SIM) in order 
to generate the data vector for each point. The data vector, which is partly the generated sig- 
nal vector and partly a vector of "ground truth” proportions, is packed by a special subroutine. 
When a complete packed line of data has been generated, it is passed back to OUT94 to be 
written out on tape. 

In order to interface with Assembly Line, External Function SIM contains the Assembly 
Line Common Block of Q-variables. In the first part of the program certain of these Q-vari- 
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ables are set, and certain SIM variables are set to zero or their default values. The variables 
SCAL, MODE, IPRINT, SUNIT, IPR1N2 are then read via a READ AND PRINT DATA statement, 
and a data title is read in. The signatures for the '’user” and "alien" materials are read in via 
READMX, with SCAL = 51,2 if CDC-1604 signatures are used. Finally, the are computed — 
the first six for "user" signatures and the last three for the "alien" signatures. These are 
normalized to sum to 1. The arguments for the random number generators (RND for uniform 
and RNP for normal distribution) are initialized to zero and SIM then returns to the enabling 
program. 

The remaining portion of SIM consists entirely of Internal Function LINE, which generates 
a complete packed line of multispectral data. This is a program entry to SIM and is called by 
OUT94 each time a line is to be generated. If the line number requested is too large, LINE re- 
turns to OUT94 where an end-of-file is written on the tape. Otherwise it starts by setting up 
the arguments required for PAK so that the packed output line can be stored in lower Erasable. 
The random number generator is "warmed up" by exercising it 50 times. LINE then enters a 
loop which generates the individual data points for the line requested by OUT 94. This loop, 
which ends at ENDL, does the following: 

(1) Determines proportion £ of alien material in the pixel. 

(2) Selects number of "user" materials that will occur in the pixel. 

(3) Selects number of "alien" materials that will occur in the pixel. 

(4) Randomly chooses the subset of materials for (2). 

(5) Randomly chooses the subset of materials for (3). 

(6) Generates a partial proportion vector, X y , corresponding to "user" materials. 

(7) Generates a partial proportion vector, X A , corresponding to "alien" materials. 

(8) From (1), (6), and (7), computes the total proportion vector, X = (1 - |)\ u + & ) A - 

(9) Using X and the signature statistics A. and M., computes the mixture signature mean 
A and covariance M. 

(10) Chooses a sample X from the normal distribution determined by A and M. (This 
requires a Cholesky decomposition and matrix multiplication.) 

(11) Forms combined vector (x, X) and "packs" it in ERIM output format using PAK. 

(12) Goes to next point; if last point is done, calls final entry to PAK. 

(13) Returns to OUT94 and writes packed line on tape. 

The number of components for X is always 9, although some may be zero. The number of data 
channels may be as large as 15, so the total vector of data plus ground truth proportions may 
have as many as 24 components on the output. When OUT94 has written the last line, it returns 
to the enabling program, which then terminates. 
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How to Use SIM 

The program deck for generating simulated data must include the following: 

• A source deck for the enabling program (described below) 

• Binary decks for SIM, PAKH, and READMX 

• Data (described below) 

The enabling program must be of the form 

*MAIN PROGRAM TO ENABLE SIM 
N'R 

P'N Q(225) 

"Dimension," "Erasable" statements 
F’T__, 

READ AND PRINT DATA 
Mount signature tape on SUNIT 

SIM.(NSIG, NCHAN, NSS, TAU1, TAU2, NUSER, NALIEN, ALFA, BETA, GAMMA) 

OUT94. (F, L, S, UNIT, BIN1, BIN2) 

E’M 

and the arguments to SIM and OUT94 are read in via the READ AND PRINT DATA statement. 
SIM itself has certain "input" variables which need not be read in since they are set to default 
values. These are 

SCAL (default = 1.0) 

MODE (default = 1) 

IPRINT (default = OB 
IPRIN2 (default = OB 

and they can be changed via the READ AND PRINT DATA statement of SIM. 

A list of input variable definitions follows: 

NSIG = number of signatures used (NALIEN and NUSER) 

NCHAN = number of data channels to be used 
NSS = number of points per line generated by SIM 
TAU1 = value of rfor computing pj lor "user” signatures 
TAU2 = value of rfor computing p. for "alien" signatures 
NUSER = number of user signatures to be read 
NALIEN = number of alien signatures to be read 

ALFA, BETA, GAMMA = parameters used to calculate £ (defined previously) 

F = first line number 
L = last line number 
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S - line number interval 

SUNIT = tape drive for signature (input) tape 

UNIT = tape drive for output tape 

BIN1, BIN2, bin location of first and last output tapes 

SCAL = scale factor for reading signatures via READ MX (should be 51.2 for 1604 tape) 
MODE = MODE argument of READMX 

IPRINT, IPRINT2 = Boolean switches for providing debugging output from SIM 

The data deck consists of a $DATA card plus input corresponding to two READ AND PRINT 
DATA statements. The first is from the enabling program and specifies values for the vari- 
ables 

NSIG, NCHAN, NSS, F, L, S, UNIT, BIN1, BIN2, TAU1, TAU2, NUSER, NALIEN, ALFA, 
BETA, GAMMA 

The second is from SIM itself and may specify any of the following: 

SCAL, MODE, IPRINT, IPRIN2, SUNIT 

which will be set to default values if not specified in the output. After these, the input deck must 
include a card specifying a title for the simulated data. The title may occupy up to 72 columns. 
An output tape with a write ring inserted into the back of it must be mounted on the tape drive 
whose number is specified by the variable UNIT. Data is generated by SIM at a rate of about 
0.5 sec per point. 

C.3 GENERAL REMARKS 

This program can provide data for testing any of several recognition or estimation tech- 
niques, even though it does notprovidefor ageometric configuration "fields.” It might be de- 
sirable to modify SIM so that more than one data point can be sampled from each mixture dis- 
tribution. 
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